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Abstract — Field  computation  by  the  stochastic  differential 
equation  (SDE)  method  is  demonstrated  for  electrostatic 
and  electrodynamic  propagation  problems  by  consider¬ 
ing  simple  examples.  The  solution  to  the  inhomogeneous 
Helmholtz  equation  is  first  related  to  that  a  Schrttdinger 
type  of  equation  (parabolic  in  nature)  by  means  of  Laplace 
transformation.  The  SDE  method  is  directly  applied  to  this 
parabolic  equation.  Presence  of  the  imaginary  term  in  the 
parabolic  equation  warrants  analytic  continuation  into  the 
complex  space  that  is  addressed  in  this  paper. 

I.  Introduction 

SDE  methods  are  alternative  methods  that  are  available 
for  field  computation  and  have  not  been  explored  to  the 
full  extent  yet.  Initial  efforts  for  solving  the  Helmholtz 
equation  are  given  in  [1]  and  [2],  although  the  approach 
of  the  latter  is  limited  to  low  wavenumbers  and  that  of 
the  former  is  geared  towards  determining  the  transport 
amplitude,  having  found  the  eikonal  by  some  other 
means.  Among  the  principal  advantages  of  the  stochastic 
method  considered  in  this  paper  are  that  (i)  it  requires  no 
meshing,  (ii)  the  field  at  a  point  can  be  determined  with¬ 
out  the  knowledge  of  the  field  at  the  neighboring  points, 
(iii)  it  is  ideal  for  parallel  computation,  and  (iv)  it  is  very 
stable  for  low  frequencies.  The  primary  disadvantage  of 
the  method  is  that  it  is  not  computationally  effective  on 
serial  machines  compared  to  the  traditional  methods.  We 
first  discuss  a  brief  theory  of  the  method  and  demonstrate 
it  for  electrodynamic  case  by  considering  the  SchrOdinger 
type  of  parabolic  wave  equation  that  is  encountered  in 
wave  propagation  problems  over  terrain. 

II.  Theory 


A.  Brownian  Motion 

A  one-dimensional  normalized  Brownian  motion  or 
Wiener  process  is  a  continuous  stochastic  process  Wt 
(W  as  a  function  of  t),  t  >  0  describing  the  motion 
of  a  particle  in  a  dynamical  system  and  satisfying  (i) 


Wo  =  0,  (ii)  for  any  0  <  t0  <h  <  •  •  •  <  tn,  the  random 
variables  Wtk  —  Wtk_1  are  independent,  1  <  k  <  n,  and 
(iii)  if  0  <  s  <  t,  Wt  —  Ws  is  normally  distributed  with 
£[Wt  -  Ws]  =  0,  £[Wt  -  Ws}2  =  (t-  s),  where  £ 
stands  for  expectation  [5].  Equivalently,  the  transitional 
probability  density  function  for  a  particle  starting  at 
position  x  and  ending  up  at  position  y  in  a  time  interval 
t  is  p{t;x,y)  =  ,x,y  6  R\t  >  0 

and  \x  —  y\  is  the  Euclidean  distance  between  x  and  y. 
Hence  the  variance  of  the  particle  displacement  increases 
linearly  with  time.  An  r-dimensional  Brownian  mo¬ 
tion  Wt  =  [W(, . . .,  Wf](  where  '  denotes  transpose, 
consisting  of  r  independent,  one-dimensional  Brownian 
motions  is  defined  with  respect  to  the  time-homogeneous 
transitional  density 

1  \x-y\2 

p{t;  x,  y)  =  /2  e  «  ,  x,  y  €  Rr,t  >  0,  (1) 


B.  The  Ito  Formula  and  the  Feynman-Kac  Formula 


If  Xt  =  satisfies  the  SDE 

dXt  =  Xt)  dWt  4-  xl)(t,Xt)dt,  where  the  ma¬ 
trix  3>(£,x)  —  {®i(tix)}  and  the  vector  = 

(t ,  x)},  then  for  a  function  f(t ,  x)  differentiable  once 
in  t  and  twice  in  x,  Ito  formula  states  that  [3] 


df{t,Xt)  = 


§-t  +  l  E 

j,k,l= 1 


_f  -  +  y  ^JL 

dxidxk  ^  dxi 
3= i 


f(t,Xt)dt 


+  E 

j,k— 1 


where  we  have  used  the  following  results:  £{dWl)  = 
0,  £(dWtjdt)  =  0,  £(dWtj)2  =  dt,  £{dWldwf)  = 
0 ,j  i=-  k.  Ito  formula  applied  directly  to  the  solution 
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u(t,x)  of  the  initial  value  problem 

du  1  r~^  a  f  .  <92  U  ^  ii/  \  &U 

at  =  + 


j,fc=i 

+  c(cc)u, 


<9a^ 


j=l 

«(0,*)  =  /(*),  (3) 

where  the  matrix  {o|.}  =  together  with  the  SDE 
dXf  =  $(Xf )dWs  +  *l>(Xf)dt,  Xjf  =  x  G  Rr 
leads  to  the  Feynman-Kac  formula 

(  t 


u(t,  a?)  —  Sx 


f(xf)  exp  |  J  c(X  f) 


ds 


(4) 


where  Sx  is  the  expectation  operator  conditioned  on 
keeping  x  fixed. 


C.  Initial  Boundary  Value  Problem 


If  the  problem  in  (3)  is  supplemented  with  the  Dirichlet 
condition  u{t ,  x)  =  g(t ,  x)  on  the  boundary  dD  of  a 
closed  domain  D ,  then  application  of  the  Feynman-Kac 
formula  to  the  extended  space-time  boundary  t  =  0  plus 
dD  yields 


u(£,  x'j  —  Sx 


Zn=t  f(Xl’x)  exp 


+^n^t  g{TUX ,  x%)  exp 


(5) 


where  E[.j  is  the  indicator  function  for  the  set  [•],  rt,x  = 
min{s  :  X*’21  G  dD },  r*  =  min(f,  Tt,x),  Xq,x  —  x , 
and  the  time  counter  for  the  stochastic  process  is  taken 
as  t\'x  —  t  —  s.  The  conditions  under  which  formulas 
(4)  and  (5)  hold  are  discussed  in  detail  in  [3].  It  is  also 
reasonable  to  expect  these  formulas  to  be  valid  for  com¬ 
plex  spatial  coordinates  assuming  analyticity  of  fields 
[1].  Essentially,  what  (5)  says  is  that  one  first  starts  a 
random  process  X\x  at  (t ,  x)  having  a  drift  component 
dictated  by  i/>  and  a  Brownian  component  dictated  by 
the  matrix  \F.  If  the  process  hits  the  boundary  dD  at 
rt,x  before  the  time  t  at  which  the  solution  in  sought, 
a  partial  contribution  comes  from  the  boundary  value 
X^xx).  On  the  other  hand,  if  the  boundary  is 
not  intercepted  before  t ,  then  a  partial  contribution  comes 
from  the  initial  value  f{X\,x).  The  process  is  then 
repeated  for  multiple  realizations  starting  at  (t ,  x)  and  an 
average  taken  to  yield  the  solution  u(t ,  x).  Even  though 
the  governing  equations  encountered  in  electrostatics  and 
electrodynamics  do  not  necessarily  take  the  form  of  (3), 
they  can  be  made  to  resemble  it  [4]  by  employing  Laplace 


transformation.  For  example,  the  electrostatic  potential 
u(x)  that  satisfies  the  Poisson  equation  V2u  =  —f(x) 
together  with  the  boundary  value  x )  can  be  written 
in  terms  of  an  auxiliary  function  K(t,x)  as  u(x)  = 

oo 

f  K(t ,  x)  dt.  The  equation  satisfied  by  K  for  t  >  0  and 

o 

its  initial,  boundary  values  are 


dK 

dt 


^2K,  K(t,  x  £  dD)  =  ip{x)  (6) 
K(t  =  0+,x)  =  ^f(x).  (7) 


In  this  case,  formula  (5)  yields  the  solution  of  the  Poisson 
equation  as 


(8) 


where  rx  is  the  first  exit  time  of  Xx  from  the  domain 
D  (see  Fig.  1).  Similarly,  solution  of  the  inhomogeneous 


dD 


Fig.  1.  Interior  domain  D  bounded  by  a  boundary  dD. 

Helmholtz  equation  V2u  +  k2u  =  —f  can  be  related  to 
the  propagator  K  via 

oo 

u(x)  =  J  K(t,x)eak  */2  dt,  (9) 

o 

where  the  propagator  satisfies  the  equation 

^  t  >  0,  K{  0,  x)  =  p(x)  (10) 

and  a  is  an  appropriately  chosen  complex  constant. 

III.  Example  Calculations 
A.  Electrostatic  Potential 

Figure  2  shows  the  electrostatic  potential  distribution 
inside  a  rectangular  region  where  the  potential  on  the 
left  wall  is  specified  as  10  volts  and  the  potential  at 
the  other  walls  is  specified  as  0  volts.  Comparison  is 
shown  between  the  analytical  solution  and  the  stochastic 


solution  that  is  obtained  with  Nr  =  3, 000  realizations. 
The  random  process  is  approximated  as  X^+At  =  Xf  + 
y/K 1 7,  where  7  is  a  zero  mean,  unit  variance  Gaussian 
random  variable.  A  good  agreement  is  seen  between  the 
two.  The  electrostatic  solution  is  rather  insensitive  to  the 
time  step  At  and  the  solution  shown  uses  At  =  0.25 
units.1 


(a)  Analytical  Solution 


Nr  =  3,000,  At'2  =  min  [Ax(r),  Ay(r),  0.5] 
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(b)  Stochastic  Solution 

Fig.  2.  Potential  distribution  inside  a  square  region. 
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B.  Propagation  Inside  a  Parallel  Plate  Waveguide 


For  the  electrodynamic  case,  we  will  consider  a  problem 
where  an  equation  of  the  type  (10)  occurs  naturally 
when  (i)  paraxial  propagation  is  considered,  and  (ii) 
back  scattering  is  neglected.  Consider  propagation  of 
time-harmonic  waves  inside  an  infinite  parallel  plate 
waveguide  with  perfectly  conducting  boundaries  at  x  —  0 
and  x  =  L.  This  model  problem  approximately  describes 
propagation  over  a  conducting  ground  with  a  periodic 
condition  enforced  at  a  large  height  L.  The  more  general 
case  of  an  impedance  boundary  condition  [6]  can  be  han¬ 
dled  using  the  same  ideas.  Let  the  axis  of  the  waveguide 
be  designated  as  t.  Assuming  an  elUT  convention  in  the 
time  variable  r  and  paraxial  propagation,  the  reduced 
field  u(t ,  x)  =  v(t ,  x)  exp  (— ikot)  satisfies  the  standard 
parabolic  equation 

du  i  d2u 
dt  2kodx2’ 

where  v(t ,  x)  is  the  true  field  and  k0  is  the  wavenumber. 
If  u(0,x)  =  uq(x),  then  the  field  at  any  distance  along 
the  guide  with  kn  =  nir/L  is  given  by 

oo 

5>(fc„)e-^t/2fc0  sin(knx),  (12) 

n— 0 

L 

2  f 

—  /  uo{x)  sin(knx)  dx.  (13) 

o 

If  the  upper  limit  in  the  summation  in  (12)  is  truncated 
to  N  and  the  integral  is  approximated  by  a  Riemann  sum 
on  a  uniform  grid  of  size  Ax  =  L/N,  then  the  relation 
between  uo(nAx)  and  Uo(kn)  will  be  that  of  Fourier 
sine  series  [6].  The  expression  for  u(t,x )  in  (12)  can 
be  analytically  continued  to  complex  space  by  writing 
x  —  £  +  ip.  The  boundaries  x  =  0  and  x  =  L  can 
also  be  analytically  continued  to  x  =  0  +  irj,  x  =  L  + 
irj  respectively.  The  field  on  the  analytically  continued 
upper  and  lower  boundaries  is  then 

DC 

u(t,  0  +  iri)  =i^2  U0{kn)e-lk*t/2ko  sinh(fc„r/)  (14) 

n=0 


u{t,  X )  = 
U0(kn)  = 


u(t,  L  +  irj)  =  i  YJ{-V)nU0(kn)e-<t/2ko  sinh 

n— 0 

(15) 

which  are  seen  to  depend  on  the  axial  distance  t.  Also 
note  that  the  field  is  no  longer  the  same  on  the  upper  and 
lower  boundaries  in  the  complex  plane.  For  the  current 
problem  it  is  desirable  to  find  the  solution  u(T,  £)  (see 
Fig.  3).  The  random  process  is  started  at  (£,T)  and  is 


Fig.  3.  Analytical  continuation  of  parallel  plate  geometry  with  three 
typical  trajectories. 

observed  over  a  duration  T.  Equation  (11)  implies  the 
governing  equation  for  the  random  process  as  dXj ^  = 
y/i/kodWs,  with  the  initial  condition  Xq  ^  =  £.  Fig¬ 
ure  3  shows  three  sample  paths  of  the  random  process. 
The  red  trace  does  not  hit  any  boundaries,  while  the  green 
and  blue  traces  hit  the  upper  and  lower  boundary  in  the 
complex  plane.  If  f(x0)  =  u( 0,£o  +  ir)0),  Pi  (T,  771)  - 
U(T  —  7i,  0  +  irfi),  and  #2(^2)  =  u(T  -  r2,  L  +  ir]2), 
then  the  solution  can  be  expressed  as 

u(T,x)  =  £e|^s<ET  f{xo)  +  ^sGti  gi{T,  T]i) 

+  Sser2  92(T,  7/2)|  (16) 

Fig.  4  shows  the  real  and  imaginary  parts  of  the  analytical 
solution  ue  and  the  stochastic  solution,  un  obtained  using 
(16)  for  a  Gaussian  initial  field  with  peak  at  x  =  Ht 
and  standard  deviation  ax.  Other  parameters  used  in  the 
computation  are  shown  in  the  figure  inset.  The  accuracy 
near  £  =  0  and  £  =  L  can  be  improved  by  optimizing  the 
time  step  At  and  the  number  of  realizations  Nr.  In  this 
example,  we  used  the  exact  solution  in  (12)  to  provide 
us  with  the  boundary  conditions  on  the  analytically 
continued  geometry  in  the  complex  plane,  and  as  such, 
the  SDE  solution  is  redundant.  However,  we  use  this 
approach  solely  to  demonstrate  the  ideas  of  the  SDE 
method.  In  more  general  case,  other  means  may  have 
to  be  found  for  specifying  the  boundary  conditions  on 
the  analytically  continued  boundaries. 

IV.  Summary 

A  theory  of  solving  the  Laplace  and  Helmholtz  equations 
using  the  SDE  approach  has  been  presented.  Favorable 


Fig.  4.  Field  at  a  distance  of  50A  from  an  aperture  source. 

comparisons  for  potential  calculations  in  electrostatics 
and  field  calculations  in  electrodynamics  have  been 
shown  by  analytical  continuation  methods.  The  electro¬ 
static  solution  is  rather  insensitive  to  the  time-stepping 
At  used  for  Brownian  motion,  but  the  electrodynamic 
case  needs  more  careful  choice.  More  complex  problems 
will  be  considered  in  the  future. 
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Introduction 


The  parabolic  equation  has  been  shown  to  accurately  model  electromagnetic  fields 
in  tunnels  for  waves  which  travel  within  ±15°  to  the  axis  of  propagation  [1].  The 
Crank-Nicolson  method  is  a  popular  finite  difference  method  that  has  been  used  to 
numerically  model  radio  wave  propagation  in  tunnels.  The  popularity  of  the  CN 
method  is  due  to  the  fact  that  it  is  stable  for  any  discretization  in  the  transverse 
plane  or  along  the  propagation  axis  [2].  The  major  limitation  of  the  CN  method  is 
that  it  requires  the  solution  of  sets  of  simultaneous  equations  that  may  become  too 
large  to  efficiently  solve  for  problems  with  dense  meshes.  The  alternate  direction 
implicit  (ADI)  technique  is  another  unconditionally  stable  FDM  which  addresses 
the  problem  of  computational  efficiency.  This  paper  presents  results  showing  the 
accuracy  of  the  ADI  technique  when  used  to  model  the  parabolic  equation  for  (a) 
square  ( b )  circular  and  (c)  semi-circular  cylindrical  PEC  tunnels.  For  each  tunnel, 
we  compare  the  numerical  solution  with  the  known  anlytical  solution  for  different 
descretizations  along  the  transverse  plane  and  propagation  axis. 

Alternate  Direction  Implicit  Method 


The  alternate  direction  implicit  method  is  a  modification  of  the  Crank-Nicolson 
method  that  can  be  used  to  numerically  solve  for  the  scalar  parabolic  equation  [2] 

du  1  /  d2  d2  \ 

dz  2jk0  \cU2  dy2  ) 

where  u(x,  y ,  z)  is  the  reduced  plane  wave  solution  and  is  related  to  the  scalar 
potential,  given  by  ^(x^y^z)  —  u(x,  y,  z)e~^koZ  [3]. The  ADI  method  reduces  the 
discretized  two  dimensional  problem  into  a  succession  of  many  one  dimensional 
problems  [2]  with  the  formulation 


un+1/2 


(2) 


4  jkQy 


n+1  _ 


=  i  + 


4  jkQ 


u 


n+l/2 


(3) 


where  (SxUm:l  ^m+l,Z  2 Umi  T  ^m— l,/)j  {.fiy'U'mJ  1  2 Umi  T  ^m,l— l)  <md 

un  represents  the  known  field,  un+1/2  is  the  unknown  virtual  field,  and  un+l  is  the 
unknown  physical  field  located  at  mAx  and  l Ay  in  a  cartesian  mesh.  Combined, 
Equations  (2)  and  (3)  are  known  as  the  Peaceman-Rachford  equations.  The  un¬ 
knowns  of  the  intermediate  plane  can  be  solved  using  Ny-1  matrices  of  rank  (A^  +  l), 
and  the  unknowns  of  the  n  +  1th  plane  can  be  solved  using  Nx-1  matrices  of  rank 
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(. Ny  +  1).  Using  the  Crank-Nicolson  method,  the  size  of  the  matrix  generated  will 
be  (Nx  +  1  )(Ny  +  1)  and  its  elements  will  not  be  in  tridiagonal  or  in  banded  form. 
As  a  result,  the  ADI  technique  will  be  more  efficient  and  solve  dense  mesh  problems 
faster.  One  difficulty  with  the  ADI  technique  is  that  the  virtual  un+l /2  field  has 
special  boundary  conditions  which  may  not  be  the  same  as  the  physical  boundary 
conditions.  For  our  following  examples,  however,  we  employ  the  same  boundaries 
for  the  virtual  and  physical  fields. 

Waveguide  examples  using  the  ADI  method 

In  the  following  sections,  we  will  use  the  ADI  method  to  solve  for  fields  in  square, 
circular  and  semi-circular  cylindrical  waveguide  type  tunnels  with  Dirichlet  and 
Neumann  boundary  conditions.  The  waveguides  operate  at  a  frequency  of  3 GHz 
and  have  a  unit  strength  gaussian  initial  field  placed  at  the  center  with  a  standard 
deviation  of  3.5A.  The  square  waveguides  have  40A  x  40A(4  x  4m)  cross  sections 
and  the  radius  of  the  circular  waveguides  are  20A(2m).  The  distance  of  propagation 
for  each  example  is  lOOOA(lOOm).  Due  to  the  Nyquist  sampling  theorem,  and  the 
limitation  of  the  PE,  our  resolution  must  be  bounded  by  Ax  —  Ay  —  2 A  [4].  The 
figures  show  pseudocolor  plots  of  the  magnitudes  of  the  analytical  and  numerical 
scalar  potentials  for  waveguides  with  Dirichlet  and  Neumann  boundary  conditions. 
The  spacings  along  the  transverse  plane  are  Ax  =  Ay  =  0.4 A  and  the  propagation 
step  size  is  Az  =  10A.  Table  1  displays  the  rms  error,  defined 


where  u^M  is  the  approximated  discretized  field,  is  the  known  analytical 

field  and  N  is  the  total  number  of  unknowns.  Figures  1  and  2  show  good  agreement 
between  the  analytical  and  numerical  fields  for  the  rectangular  and  circular  tunnels 
with  both  Dirichlet  and  Neumann  boundary  conditions.  For  the  rectangular  tunnel, 
the  Neumann  boundary  conditions  were  approximated  using  2nd  order  accurate  one¬ 
sided  approximations  and  for  the  circular  tunnel,  the  Neumann  boundary  conditions 
were  approximated  by  first  order  interior  interpolation.  In  Figure  3,  the  Neumann 
boundary  condition  were  approximated  with  first  order  one  sided  approximations. 
As  with  the  previous  two  cases,  the  field  patterns  are  closely  matched,  however  for 
the  Neumann  case  the  rms  error  is  larger. 
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Conclusions 

The  ADI  technique  for  the  electromagnetic  propagation  in  PEC  waveguides  for 
Dirichlet  and  Neumann  boundary  conditions  is  studied  and  it  is  found  that  the  field 
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patterns  of  the  numerical  solutions  are  closely  matched  with  that  of  the  analytical 
solutions.  Future  work  will  involve  the  more  irregular  shapes  and  mixed  boundary 
conditions. 
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Table  1:  The  normalized  rms  error. 


Erms  (%) 

Rectangular  WG 

Circular  WG 

Semi-circular  WG 

Ax,  Ay 

Az 

Dirichlet 

Neumann 

Dir. 

Neu. 
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Neu. 

0.8A 
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Figure  1:  The  (a)  Analytical  solution  and  (b)  Numerical  approximation  of  the  rect¬ 
angular  waveguide  with  Dirichlet  boundary  conditions  (RMS  error  =  7.3%)  and  the 
(c)  Analytical  solution  and  (d)  Numerical  approximation  of  the  rectangular  waveg¬ 
uide  with  Neumann  boundary  conditions  (RMS  error  =  7.1%)  . 
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Figure  2:  The  (a)  Analytical  solution  and  (b)  Numerical  approximation  of  the  cir¬ 
cular  waveguide  with  Dirichlet  boundary  conditions  (RMS  error  =  7.0%)  and  the 
(c)  Analytical  solution  and  (d)  Numerical  approximation  of  the  circular  waveguide 
with  Neumann  boundary  conditions  (RMS  error  =  11.9%). 


Figure  3:  The  (a)  Analytical  solution  and  (b)  Numerical  approximation  of  the  semi¬ 
circular  waveguide  with  Dirichlet  boundary  conditions  (RMS  error  =  7.6%)  and 
the  (c)  Analytical  solution  and  (d)  Numerical  approximation  of  the  semi-circular 
waveguide  with  Neumann  boundary  conditions  (RMS  error  =  27.5%). 
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Propagation  Prediction  in  Rough  and  Branched 
Tunnels  by  the  ADI-PE  Technique 


R.  Martelly* 


R.  Janaswamyt 


Abstract  —  The  study  of  radiowave  propagation  in 
tunnels  is  important  for  the  development  of  telecom¬ 
munication  systems.  The  vector  Parabolic  Equa¬ 
tion  (PE)  and  the  alternate  direction  implicit  (ADI) 
technique  are  used  to  study  radiowave  propagation 
in  branched  tunnels  and  in  tunnels  with  rough  walls. 
Previous  work  has  shown  that  the  ADI-PE  method 
can  accurately  predict  transmission  loss  in  straight 
tunnels  with  smooth  walls  and  with  known  electrical 
parameters.  We  extend  the  analysis  of  this  method 
by  including  the  realistic  cases  of  branching  tunnels 
and  tunnels  with  rough  walls.  We  breifly  discuss 
the  boundary  conditions  used  in  each  case  as  well  as 
compare  our  results  with  known  numerical  or  ana¬ 
lytical  models.  The  numerical  results  obtained  for 
the  branched  tunnel  case  were  compared  with  the 
results  produced  by  HFSS  and  the  results  for  the 
rough  wall  case  were  compared  with  known  analyt¬ 
ical  loss  factors.  The  results  show  excellent  agree¬ 
ment  in  both  cases. 


than  the  Crank-Nicolson  method  [1].  ADI  intro¬ 
duces  slightly  more  error  than  the  Crank-Nicolson 
method  but  previous  work  has  shown  that,  for  mod¬ 
est  discretizations,  the  ADI  and  Crank-Nicolson  so¬ 
lutions  are  nearly  identical  [1].  After  a  brief  discus¬ 
sion  on  the  ADI-PE  formulization  (section  2),  we 
continue  our  analysis  of  radiowave  propagation  in 
tunnels  using  the  ADI-PE  by  considering  the  spe¬ 
cial  cases  of  branched  tunnels  (section  3)  and  tun¬ 
nels  with  rough  walls  (section  4). 

2  ADI-PE  Theory 

The  vector  PE,  as  formulated  by  Popov  [2],  defines 
the  transverse  electric  fields  in  terms  of  a  vector 
function,  W,  as  shown  in  equation  (1)  [2] 


1  INTRODUCTION 


(Ex,Ey)T  =  We~ik°z,  (1) 


The  alternate  direction  implicit  (ADI)  method  cou¬ 
pled  with  the  vector  parabolic  equation  (PE)  has 
previously  been  shown  to  model  radiowave  propa¬ 
gation  in  straight  tunnels  with  smooth  walls  [1,  2]. 
The  PE  model  assumes  the  propagating  fields  are 
slowly  varying  in  the  direction  of  propagation  and 
backscattered  fields  are  ignored.  In  realistic  tun¬ 
nels,  over  a  long  distance,  higher  order  propagating 
modes  attenuate  and  the  lowest  order,  slowly  vary¬ 
ing  modes,  dominate  [3].  The  vector  PE  can  ac¬ 
curately  solve  for  the  lower  order  dominate  modes 
which  travel  within  ±15°  to  the  axis  of  propaga¬ 
tion. 

The  vector  PE,  as  formulated  by  Popov  [2],  char¬ 
acterizes  the  electrical  parameters  of  the  tunnel 
walls  with  an  equivalent  surface  impedance  and  en¬ 
forces  the  impedance  boundary  condition  [2].  The 
Crank-Nicolson  method  has  been  traditionally  used 
to  numerically  solve  the  vector  PE  because  it  is 
an  unconditionally  stable  FDM;  but  it  can  also  be 
computationally  intensive  [1].  The  ADI  technique 
is  also  an  unconditionally  stable  implicit  FDM 
that  is  significantly  more  computationally  efficient 
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where  kQ  is  the  free  space  wave  number  and  z  is  the 
direction  of  propagation.  The  vector  function,  W, 
describes  the  complex  wave  amplitude  of  the  plane 
wave.  The  vector  PE  can  then  be  defined  as 


2  Kj 


W 

~dz 


d2W  d2W 
dx 2  dy 2 


(2) 


When  the  impedance  boundary  condition  is  en¬ 
forced  on  the  tunnel  wall,  the  transverse  E  fields 
are  coupled  and  the  value  of  W  at  the  boundary  is 
shown  to  be 


w=AJn*  nv  Y l!Zs  0 

k\  ^ y  ^ x 


nx 

fly 


ny  \dW 
dfi 


(3) 


where  nx  and  ny  are  the  x  and  y  components  of 
the  normal  vector  on  the  boundary  and  Zs  is  the 
normalized  surface  impedance. 

As  shown  in  [1],  the  ADI-PE  can  be  summarized 
by 


where  r^)  is  the  mesh  ratio  and  the  difference  quo¬ 
tient,  Sx(y),  represent  the  second  order  discretiza¬ 
tion  of  the  x  and  y  derivatives. 
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3  Branched  Tunnels 

3.1  Branched  Tunnel  Model 

A  typical  branch  tunnel  geometry  is  shown  in  Fig¬ 
ure  1,  where  the  main  tunnel  axis  is  shown  as  a 
solid  bold  line  and  the  branch  tunnel  axis  is  shown 
as  the  long  dashed  line.  The  branch  angle,  0*,,  is  the 
angle  between  the  axis  of  the  straight  and  branch 
tunnel.  Using  this  geometry,  the  branch  angle  must 
be  less  than  30°  to  ensure  that  reflected  rays  enter¬ 
ing  the  branch  are  within  our  PE  limit  of  ±15°. 
Considering  also  the  rays  diffracted  from  the  cor¬ 
ners  at  the  junction,  we  can  arrive  at  a  much  more 
stringent  requirement  of  Ob  <  15.  Figures  1  and  2 
show  the  incident  and  reflected  ray  as  it  enters  the 
branch  tunnel.  The  short  dashed  line  markes  the 
input  plane  of  the  branch  tunnel.  The  grazing  an¬ 
gle,  T,  and  the  angle  of  the  ray  entering  the  branch, 
a ,  are  also  shown. 


Figure  1:  The  incident  and  reflected  ray  entering 
the  branch  tunnel  when  T  <  0&. 


Figure  2:  The  incident  and  reflected  ray  entering 
the  branch  tunnel  when  T  <  0^. 

The  slope  of  the  branching  wall  is  modeled  using 
a  staircase  approximation  (see  Figure  3a)  and  the 
impedance  boundary  condition  is  enforced  on  all 
four  walls  as  outlined  in  [1]  and  [2].  The  fields  along 
the  planes  marking  the  entrance  of  the  main  tunnel 


(line  C)  and  the  branch  tunnel  (line  B)  are  solved 
simultaneously  and  then  used  as  the  initial  fields  for 
the  two  separate  diverging  tunnels.  We  simulate  a 
0.9  mx  1.0  m  rectangular  tunnel  with  a  branch  angle 
of  15°  and  operating  at  a  frequency  of  900  MHz.  We 
used  the  far  zone  field  of  a  unit  strength  Gaussian 
field  source  as  our  initial  field.  The  source,  with 
standard  deviation  of  0.75  A  (wavelength),  is  placed 
at  the  center  of  the  tunnel  entrance.  This  means  we 
are  only  using  the  lowest  order  modes  as  our  initial 
field  [1,  3].  The  fundamental  mode  propagates  near 
our  PE  limit  at  14.17°. 

3.2  Comparison  of  numerical  results  to 
HFSS 

For  our  ADI  simulation,  we  used  discretizations 
(within  the  tunnel  junction)  of  Ax  =  0.060  A, Ay 
=  0.054  A  and  Az  —  0.013  A.  The  cross-sectional 
coordinates  are  indicated  by  x  and  y,  while  the  ax¬ 
ial  coordinate  in  the  main  tunnel  is  denoted  by  the 
2:- axis.  To  validate  our  results,  we  compared  our  so¬ 
lutions  with  HFSS  [6]  and  plot  the  magnitude  of  the 
Ey  field  along  the  main  and  branch  tunnel  axis.  In 
the  HFSS  simulation,  we  used  radiation  boundary 
conditions  to  terminate  the  tunnel  and  symmetry 
planes  to  reduce  computational  labor. 

HFSS  is  a  full  wave  simulator  and,  unlike  the 
ADI-PE,  solves  for  backscattered  waves  as  well 
as  waves  traveling  in  the  forward  direction.  The 
backscattering  is  seen  as  fluctuations  in  the  axial 
field  in  Figure  3b  near  the  diverging  tunnels.  The 
tunnel  has  a  dielectric  constant  of  5  and  conduc¬ 
tivity  of  0.1  S/m.  A  high  conductivity  was  chosen 
so  there  will  be  appreciable  loss  in  the  small  tunnel 
dimensions  allowed  in  HFSS.  As  we  can  see  from 
Figure  3b,  there  is  strong  agreement  in  the  axial 
field  intensity  along  the  main  and  branch  tunnel 
axis  between  the  ADI  and  HFSS.  The  figure  also 
shows  that  there  is  about  a  10  dB  drop  when  going 
from  the  main  to  the  branch  tunnel  (at  the  point 
marked  C  in  Figure  3b). 

Although  the  ADI-PE  was  used  to  simulate 
a  tunnel  with  a  relatively  small  electrical  cross- 
section  (2.7Ax3A),  the  ADI-PE  is  capable  of  han¬ 
dling  larger  tunnels  at  higher  frequencies  without 
running  into  memory  problems  on  an  average  (2  GB 
RAM)  PC  [1]. 

4  Tunnels  with  Rough  walls 

4.1  Surface  roughness  model 

Surface  roughness  is  the  local  variation  of  the  tun¬ 
nel  wall  relative  to  a  mean  surface  level  [4,  5].  In 
this  study  we  consider  random  surface  deviations 
in  an  otherwise  smooth  wall.  For  the  purpose  of 
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Figure  3:  (a)  Geometry  of  the  branch  tunnel,  (b)  The  axial  field  intensity  of  the  main  tunnel  from  ADI 
(blue,  solid),  HFSS  (dark  green,  dashed),  the  branch  tunnel  from  ADI  (red,*)  and  HFSS  (light  green, 
dot). 


numerical  computations  we  assume  the  random  de¬ 
viations  to  be  Gauassian  distributed.  A  Gaussian 
distribution  of  the  surface  level  can  be  character¬ 
ized  by  a  root-mean-square  height  deviation  a h  and 
correlation  length,  /  [4,  5].  Smooth  tunnels  have  a 
typical  RMS  height  deviation  of  0.01m  and  rough 
surfaces,  such  as  those  found  in  coal  mine  tunnels, 
have  a  RMS  height  deviation  of  0.1m  [4].  The  loss 
of  the  Ex  field  in  a  rectangular  tunnel  is  given  by 
equation  (5)  [4] 

Lrough  =  4.343tt2/i2A2  (T  +  ^jz  (dB)  (5) 

where  d\  and  #2  are  the  tunnel  dimensions  in  the  x 
and  y  axis,  respectively. 

Small  scale  roughness  can  be  characterized  by  re¬ 
placing  the  rough  surface  with  a  flat  impedance  sur¬ 
face  that  produces  an  equivalent  specular  reflection 
coefficient.  The  equivalent  impedance  for  horizon¬ 
tal  polarization  is  given  by  [5] 

( Zs-j(k0ah )2^,  kj«  1 

Zeq=  <  Zs-\- (k0(Jh)  sin's!/,  kQl^> l,\I/^5> 

jzs  +  (fc007j)2r(!)^/  fc0Z»l,'J/«^7= 

(6) 

where  Zs  is  the  surface  impedance  for  the  smooth 
wall,  r(-)  is  the  Gamma  function,  and  T  is  the  graz¬ 
ing  angle.  Due  to  the  PE  angle  limitation  of  ±15°, 
the  maximum  slope  angle  0  of  the  rough  surface 


and  the  grazing  angle  must  satisfy  the  following  re¬ 
lationship 

20  +  T  <  15°  (7) 

where  T  is  defined  as  shown  in  Figure  4.  As  we 
can  see  from  Figure  4,  the  angle  of  the  specular 
reflection  of  the  incident  ray,  denoted  by  £,  will 
depend  on  the  height  deviation  of  the  roughness. 
The  roughness  angle,  0 ,  is  related  to  and  l  by 
tan#  =  crfc/Z. 


Figure  4:  The  geometry  of  the  rough  surface. 

4.2  Comparison  of  numerical  and  analyti¬ 
cal  solutions 

We  consider  a  rectangular  4.26 mx2. 10m  tunnel 
and  a  circular  tunnel  with  radius  of  2  m.  The  fun¬ 
damental  EIFjq  mode  is  used  as  the  initial  field  of 
the  rectangular  tunnel  and  the  fundamental  TE0i 
mode  generated  by  a  loop  ring  excitation  is  used  as 
the  initial  field  for  the  circular  tunnel.  Both  tunnels 
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operate  at  a  frequency  of  1  GHz  and  the  dielectric 
constant  and  conductivity  of  the  tunnel  wall  are 
taken  as  er=12  and  aQ  =  0.02  S/m,  respectively. 

Tables  1  and  2  summarize  the  mode  attenua¬ 
tion  factors  (MAFs),  or  the  loss  indB/km,  of  the 
smooth  and  rough  tunnels  with  rectangular  and  cir¬ 
cular  cross-sections,  respectively.  In  real  tunnels,  as 
in  our  simulations,  the  lowest  order  mode  will  de¬ 
termine  the  MAF  over  a  long  distance.  We  used 
equation  (5)  as  our  analytical  loss  factor  for  both 
the  rectangular  and  circular  tunnel.  As  we  can  see 
from  equation  (5),  the  loss  due  to  roughness  is  a 
function  of  wavelength.  To  notice  an  appreciable 
loss  at  1  GHz,  we  need  to  assume  the  walls  are  as 
rough  as  cave  walls.  Therefore,  we  used  a  RMS 
height  of  0.1m  (0.33  A)  and  a  correlation  length  of 
2.5  m  (8.33  A)  for  both  tunnels. 

The  grazing  angle  of  the  fundamental  mode  of 
the  rectanglular  and  circular  tunnel  is  4.56°  and 
5.25°,  respectively.  The  roughness  angle  is  2.29° 
and  equation  (7)  is  satisfied.  The  ADI  is  simulated 
using  discretizations  of  Ax  =  0.284  A,  Ay  =  0.14  A 
and  Az  =  4  A  for  the  rectangular  tunnel  and  Ax  = 
Ay  =0.44  A  and  Az  =  1.67  A  for  the  circular  tunnel. 
As  we  can  see  from  Table  1,  the  excess  loss  due  to 
roughness  for  the  rectangular  tunnel  is  about  7  dB 
when  using  either  equation  (5)  or  the  ADI-PE.  Sim¬ 
ilarly,  Table  2  shows  the  same  close  agreement  in 
numerical  and  theoretical  excess  loss  for  the  tunnel 
with  circular  cross-section.  In  this  case,  the  excess 
loss  due  to  roughness  is  about  16  dB. 


Rectangular  Tunnel  (dB/km) 

4.26  mx2.10  m 

Analytical 

ADI-PE 

w/o  Roughness 

31.0 

29.5 

Roughness 

38.0 

36.6 

Table  1:  Analytical  and  Numerical  MAFs. 


Circular  Tunnel  ( dB /km) 
Radius  =  2.0  m 

Analytical 

ADI-PE 

w/o  Roughness 

11.0 

10.1 

Roughness 

27.0 

26.1 

Table  2:  Analytical  and  Numerical  MAFs. 


and  electrically  small  tunnel  cross-sections  and  the 
equivalent  Crank-Nicolson  code  would  require  sig¬ 
nificantly  larger  matrices  [1]. 

5  Conclusions 

The  ADI-PE  method  has  been  shown  to  accu¬ 
rately  model  branch  tunnels  and  tunnels  with  rough 
walls.  For  branch  tunnels,  even  at  the  PE  limit, 
there  is  good  agreement  between  the  ADI-PE  and 
with  commerical  simulation  codes  such  as  HFSS. 
Also,  the  additional  loss  created  by  tunnels  with 
rough  walls  is  correctly  modeled  using  equivalent 
impedances.  The  ADI-PE  method  compares  well 
with  known  theoretical  roughness  loss  factors. 
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The  accuracy  of  the  results  suggests  that  the 
equivalent  surface  impedance,  along  with  ADI-PE, 
is  an  adequate  model  for  determining  loss  due  to 
surface  roughness.  Unlike  the  ADI-PE,  a  finite  el¬ 
ement  method  would  be  limited  to  low  frequencies 
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Abstract 

The  diffusion  and  Schrodinger  propagators  have  been  known  to  coexist  on  a 
lattice  when  a  particle  undergoing  random  walk  is  endowed  with  two  states  of 
spin  in  addition  to  the  two  states  of  direction  in  a  1+1  spacetime  dimension. 
In  this  paper  we  derive  explicit  expressions  for  the  various  transitional 
probabilities  by  employing  generating  functions  and  transform  methods.  The 
transitional  probabilities  are  all  expressed  in  terms  of  a  one-dimensional  integral 
involving  trigonometric  functions  and/or  Chebyshev  polynomials  of  the  first 
and  second  kind  from  which  the  spacetime  continuum  limits  of  the  diffusion 
equation  and  Schrodinger  equation  follow  directly. 

PACS  numbers:  03.65.— w,  05.40.Fb 


1.  Introduction 

There  has  been  a  lot  of  interest  in  the  recent  past  to  understand  quantum  mechanics  in  the 
context  of  classical  statistical  mechanics.  On  the  one  hand,  Brownian  motion  provides  a 
microscopic  model  of  diffusion  and  provides  an  unambiguous  interpretation  of  the  diffusion 
equation.  On  the  other  hand,  a  similar  physical  interpretation  is  lacking  for  the  Schrodinger 
equation,  whose  wave  solution  is  a  complex  quantity  without  a  physical  reality.  Because 
classical  diffusion  cannot  account  for  the  self-interference  pattern  that  is  so  intrinsic  to  quantum 
behavior,  several  theories  have  been  put  forward  recently  to  understand  the  microphysics  of 
quantum  behavior.  Nelson  [1]  derived  the  Schrodinger  equation  starting  from  Newtonian 
mechanics  and  by  assuming  that  a  particle  is  subject  to  an  underlying  Brownian  motion 
described  by  a  combined  forward-in-time  and  a  backward-in-time  Wiener  processes.  A 
detailed  account  of  Nelson’s  original  idea  of  stochastic  mechanics  and  its  subsequent 
refinement  is  given  in  [2-5].  Nottale  [6]  and  Ord  [7]  advanced  the  idea  that  spacetime  is 
not  differentiable  but  is  of  a  fractal  nature,  suggesting  that  an  infinity  of  geodesics  lie  between 
any  two  points  and,  thereby,  providing  a  fundamental  and  universal  origin  for  the  double 
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Wiener  process  of  Nelson.  These  ideas  are  elaborated  in  the  monograph  [8].  El  Naschie  [9] 
too  considered  a  fractal  spacetime  with  a  Cantorial  structure  and  argued  that  quantum  behavior 
could  be  mimicked  by  combining  this  fractal  spacetime  with  a  diffusion  process.  A  totally 
different  paradigm  was  recently  introduced  by  Ord  [10],  who  by  considering  a  symmetric 
random  walk  on  a  lattice,  showed  that  both  the  diffusion  equation  and  the  Schrodinger  equation 
occur  as  approximate  descriptions  of  different  aspects  of  the  same  classical  probabilistic 
system.  By  considering  a  4-state  random  walk  (4RW)  on  a  discrete  lattice,  wherein  a  particle 
is  endowed  with  two  states  of  direction  and  two  states  of  spin,  Ord  [10-12]  has  shown  that 
both  diffusion  and  Schrodinger  propagators  coexist  on  a  lattice  and  that  either  can  be  obtained 
from  a  distinct  projection  of  the  same  random  walk.  It  is  too  early  to  speculate  as  to  which  of 
Nelson’s  or  Ord’s  model  will  duplicate  the  true  quantum  behavior  under  a  variety  of  situations. 
This  can  only  be  ascertained  through  additional  work  on  both  models.  It  may  be  mentioned 
that  the  combination  of  displacement  and  spin  have  also  been  used  previously  in  [13,  14] 
to  study  dynamics  of  a  quantum  particle  in  spacetime.  However,  the  important  distinction 
between  the  Ord  model  and  the  one  considered  in  [13,  14]  is  that  the  states  describing  the 
direction  of  motion  are  independent  of  those  describing  the  spin  states  in  the  former  model. 
There  is  also  an  intrinsic  notion  of  memory  embedded  in  the  Ord’s  model. 

The  Schrodinger  type  of  equation  is  encountered  under  the  guise  of  parabolic  wave 
equation ,  or  simply  parabolic  equation  in  the  solution  of  boundary-value  problems  in  several 
branches  of  applied  physics  such  as  acoustics  [15],  optics  and  classical  electromagnetic  wave 
propagation  [16].  In  such  boundary- value  problems,  inhomogeneities  of  the  propagating 
medium  caused  by  the  varying  index  of  refraction  of  the  intervening  material  take  the  place 
of  the  potential  field  experienced  by  a  quantum  particle.  The  standard  parabolic  equation 
is  resulted  when  one  extracts  paraxial  propagation  along  a  preferred  direction  from  the  full 
Helmholtz  equation.  In  addition  to  providing  a  microscopic  model  for  the  Schrodinger 
equation,  the  4RW  model  considered  by  Ord  is  also  attractive  in  the  solution  of  stochastic 
differential  equations  associated  with  these  parabolic  type  of  equations,  carried  out  by 
employing  only  real  random  processes.  Because  walks  modeling  the  Schrodinger  equation 
in  the  4RW  model  traverse  only  real  space,  no  analytical  continuation  of  boundary  data  into 
complex  space  is  required  that  would  otherwise  be  demanded  [17,  18]  when  solving  these 
boundary- value  problems. 

Ord  does  not  provide  explicit  expressions  for  the  various  transitional  probabilities,  but, 
instead,  discusses  the  continuum  limits  directly  from  the  governing  difference  equations.  For 
a  variety  of  reasons,  it  is  desirable  to  obtain  closed-form  expressions  (or  those  involving 
integrals)  for  these  transitional  probabilities.  In  this  paper,  we  provide  analytical  expressions 
for  the  transitional  probabilities  associated  with  the  4-state  random  walk  in  1+1  dimension 
in  spacetime  by  using  a  transform  approach.  Our  work  here  is  partly  motivated  by  the 
desire  to  have  expressions  for  the  transitional  probabilities  while  solving  the  aforementioned 
boundary- value  problems  using  the  parabolic  equation  in  a  homogeneous  medium.  Using 
these  expressions,  it  is  further  shown  that  in  the  continuum  limits  as  the  mesh  size  shrinks  to 
zero  in  both  space  and  time,  one  directly  recovers  the  diffusion  equation  and  the  Schrodinger 
equation.  Thus,  the  main  contributions  of  the  paper  are  to  (i)  elucidate  methodology  for 
obtaining  the  closed-form  expressions  for  the  various  transitional  probabilities  of  the  4RW, 
and  (ii)  establish  the  continuum  limits  of  the  diffusion  and  Schrodinger  equations  describing 
the  dynamics  of  particles  obeying  the  4RW.  The  methodology  presented  in  this  paper  is 
most  suitable  for  describing  quantum  dynamics  of  a  free-particle,  although  the  4RW  model 
itself  has  been  extended  in  the  presence  of  a  potential  field  [19].  The  paper  is  organized 
as  follows:  section  2  gives  a  brief  introduction  of  the  random  walks  considered  in  [10,  12]. 
Section  3  introduces  the  generating  functions  and  the  2D  transforms  considered  in  this  paper. 
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Table  1.  Various  states  in  random  walk. 


State 

Direction 

Spin 

1 

Right 

+1 

2 

Left 

+1 

3 

Right 

-1 

4 

Left 

-1 

Section  4  provides  expressions  for  the  various  transitional  probabilities  as  well  as  discusses 
the  derivation  of  the  diffusion  equation  and  the  Schrodinger  equation  as  continuum  limits  of 
these  probabilities. 


2.  Multistate  random  walks 

Consider  the  4RW  model  proposed  by  Ord  and  Deakin  [12],  where  a  particle  undergoes  random 
motion  in  discrete  spacetime  (x  =  mA,t  =  56),  with  x  denoting  space  and  t  denoting  time, 
and  A  and  6  denoting  the  spatial  and  temporal  steps,  respectively.  At  every  point  the  particle 
is  endowed  with  two  independent  binary  properties,  its  direction  of  motion  (right  or  left)  and 
its  spin  or  parity  (±1).  The  particle  is  assumed  to  change  its  direction  with  every  collision, 
but  change  its  spin  only  every  other  collision.  The  four  states  of  the  particle  corresponding 
to  the  four  combinations  of  direction  and  spin  are  indicated  in  table  1 .  Note  that  the  particle 
can  execute  any  direction  of  motion  irrespective  of  the  spin,  in  contrast  to  the  model  used  in 
[13,  14].  However,  there  is  an  intrinsic  assumption  of  memory  in  Ord’s  model  that  arises  from 
keeping  track  of  the  parity  of  collisions.  If  p^tm A,  56) A,  /x  =  1,  . . . ,  4,  is  the  probability 
that  a  particle  is  in  state  /x  at  the  spacetime  point  (m  A ,  56),  m  =  0,  d=l,  ±2,  . . . ,  5  =  0,  1, ... , 
then  the  transitional  relations  considered  in  [12]  were  of  the  form 

Pi[mA,  ( s  +  1)6)]  =  ap\[(m  —  1)A,  se]  +  fip\\{m  +  1)A,  56] 

P2[mA,  ( 5  +  1)6)]  =  otp2[(m  +  1)A,  56]  +  /3p\[(m  —  1)A,  56] 

(1) 

pilmA,  (s  +  1)6)]  =  apiKm  —  1)A,  56]  +  /3p2[(m  +  1)A,  56] 

P4[mA,  ( s  +  1)6)]  =  ap^iim  +  1)A,  56]  +  fipslim  —  1)A,  56], 

where  a  +  /3  =  1.  Here,  a  is  the  probability  that  a  particle  maintains  its  direction  at  the  next 
time  step,  whereas  /3  is  the  probability  that  it  will  change  its  direction  at  the  next  time  step. 
The  Markov-chain  character  of  the  transitional  probabilities  is  evident  from  definitions  in  (1). 
From  the  total  probability  theorem,  the  probability  that  a  particle  is  somewhere  on  the  lattice 
at  a  given  time  is  equal  to  1  and  is  represented  mathematically  by 

4  oo 

E  E  Pn(mA,  se)A  =  1.  (2) 

li=  1  m=—o o 

Ord  [10]  has  shown  that  the  diffusion  and  Schrodinger  propagators  coexist  on  the  lattice 
and  that  both  behaviors  are  embedded  in  equations  (1).  To  affect  a  separation  of  the 
diffusive  behavior  from  the  wave-like  behavior,  the  following  linear  transformation  is 
used:  qi(mA,se)  =  2s^2[p\(mA,  se)  —  p3(mA,s€)],q2(mA,s€)  =  2s^2[p2(mA,  se)  — 

P4(mA,  56)],  w\{m A,  56)  =  [ p\{mA ,  56)  +  p2(mA ,  56)  +  p^imA,  se)  +  p4(mA,  56)],  and 

W2(mA,se)  =  [pi(mA,se)  +  p^imA.se)]  —  [p2(mA,se)  +  p4(mA,se)].  The  quantity 
qiA  (without  the  weight  factor  2s /2)  indicates  the  expected  difference  in  the  number  of 
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particles  of  opposite  spin  arriving  at  (m A,  sc)  while  moving  to  the  right.  Similarly,  q2A 
refers  to  the  expected  number  of  particles  arriving  at  (mA,s€)  while  moving  to  the  left. 
Also,  w\{mA,s€)A  is  the  probability  that  a  particle  leaves  (mA,s€)  in  either  direction 
and  in  any  spin  state,  and  W2(mA,  se)  A  is  the  difference  in  the  probabilities  that  a  particle 
leaves  (mA,  sc)  to  the  right  and  the  left.  Introducing  the  shift  operator  Efx p^m A,  se)  — 
p^iim  =b  1)A,  s€],  a  time- advancing  operator  Etp^imA,  sE)  —  p^lmA,  (s  +  1)6],  and  the 
vector  p  =  [pi,  p2,  ps,  Pa\t  ,  where  the  superscript  T  denotes  transpose,  the  transitional 
relations  in  (1),  which  are  of  the  form  Et p  =  S*p,  get  transformed  into 


E  A1)  =  -(  {Ex  + E;V)  Ex ~ E;V)  \  A1) 

'  W  2  \(a  -  p)(Ex  -  E ;l)  (a  -  p)(Ex  +  E;1) )  W  ’ 


Thus  the  variables  (w i,  w2)  get  decoupled  from  (qi,  q2).  Essentially,  this  decoupling  results 
from  block-diagonalizing  the  matrix  S*  and  describing  the  system  in  terms  of  its  eigenstates. 
The  physical  significance  of  this  transformation  is  touched  upon  in  [11,  12].  Note  that 
wj  and  qj  need  not  strictly  be  probabilistic  quantities  (meaning  ^0),  but  we  will  continue  to 
describe  them  as  ‘transitional  probabilities’  with  the  understanding  that  the  actual  probabilistic 
quantities,  namely,  p can  be  easily  recovered  from  these  using  the  inverse  relations. 


3.  Generating  functions  and  transforms 


We  are  interested  in  the  solutions  of  (3)  and  (4)  for  the  special  case  of  a  symmetric  random 
walk  with  a  =  =  0.5.  In  this  case  we  have  a  set  of  linear  difference  equations  and 

the  solution  can  be  obtained  conveniently  using  transform  methods  [20,  21]  and  appropriate 
generating  functions.  The  key  step  here  is  to  pick  a  suitable  transform  consistent  with  the 
nature  and  domain  of  definition  of  the  problem.  We  denote  the  2D  transform  £,  consisting 
of  a  Fourier  transform  in  space  (owing  to  the  unbounded  nature  of  the  spatial  coordinate)  and 
the  z-transform  [22]  in  time  (the  z-transform  can  be  arrived  from  the  discretized  version  of  a 
Laplace  transform  and  is  suitable  for  discrete  functions  defined  on  a  half-line),  of  a  discrete 
function  v{m A,  se)  as  V ( kx ,  z)  and  define 

CO  CO 

V(kx,z)  =  A  ^  ^(^A,  se)zs  Q~lmkxA  =  Cv(mA,  se).  (5) 

m=— oo  5=0 

The  inverse  relation  can  then  be  obtained  as 

i  r*/A  r  y(b  7\ 

v(mA,se)  = J  j)  — ^ — emkxA  dkxdz  =  C~lV(kx,  z),  (6) 

where  the  identities 

f  ei(n~m)kxAdkxA  =  27tSl  (7) 

J  kxA  =—7T 


d  z  =  27ri<^ 


(8) 


are  used  to  derive  (6).  Here  5"  is  the  Kronecker’s  delta  and  Cz  is  a  closed  contour  around  the 
origin  in  the  complex  z-plane  that  encloses  only  the  singularities  at  the  origin.  The  present 
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analysis,  consisting  of  the  z-transform  along  the  time  axis  and  Fourier  transform  along  the 
spatial  axis,  is  most  suitable  for  studying  linear  difference  equations  with  constant  coefficients 
such  as  encountered  in  the  study  of  free-Schrodinger  equation  by  the  4RW  model.  Other 
suitable  methods  must  be  devised  for  studying  particle  motion  in  the  presence  of  a  potential 
field.  Note  that  V ( kx ,  z)  is  periodic  in  kx  with  a  period  2tt / A.  Using  the  definition  in  (5),  it 
can  also  be  shown  that 

Cv[mA,  (s  +  l)e]  =  z"1  [V(kx,  z)  -  V0(kx)]  (9) 

Cv[(m  ±  1)A,  se]  =  e±ik'AV(kx,  z),  (10) 

where  Vo  (kx )  is  the  Fourier  transform  of  the  initial  distribution  v(m  A,  0): 

00 

V0(kx)  =  A  J2  v(mA,0)e-“'4.  (11) 

m = — oo 

Note  that  the  periodicity  property  of  Vq (kx)  implies  that  Vq(tt/ A)  =  Vo(—tt/  A). 


4.  Transitional  probabilities 

Having  defined  the  required  transforms,  we  will  now  derive  expressions  for  the  transitional 
probabilities  uq,  u>2,  qi  and  q2.  Because  of  the  decoupling  afforded  in  (3)  and  (4),  it  is 
sufficient  to  consider  the  diffusive  and  wave-like  behaviors  separately. 


4.1.  Diffusive  behaviour 

The  diffusive  part  of  the  particle  motion  is  governed  by  the  discrete  functions  w\  and  W2  as 
will  be  evident  shortly.  Let  W\{kx,  z)  and  W2(kx,  z)  be  the  2D  transforms  of  uq(raA,  se)  and 
u;2(raA ,  56)  and  Yife)  and  T2 (kx)  be  the  transforms  of  the  initial  distributions  w\(m A,  0) 
and  iu2(raA,  0),  respectively.  From  the  definition  of  w\  in  terms  of  p^,  /x  =  1,  . . . ,  4,  and 
relation  (2),  it  is  seen  that  Yi(0)  =  1.  On  applying  the  transform  C  to  the  set  (3)  and  making 
use  of  the  properties  (9)  and  (10),  it  is  easy  to  see  that  z)  =  Y2fe)  and 


W!(kx,z) 


Ttfej-usinfeAjTzfe) 
1  —  z  cosfe  A) 


00 

=  ^  z"  cosn(kxA)  [Ti (kx)  -  iz  sinfe  A)T2fe)] 

n= 0 


(12) 

(13) 


where  (13)  has  been  obtained  by  using  the  series  expansion  of  [1  —  z  cosfe  A)]-1.  Such  a 
series  converges  uniformly  provided  that  |z  cos(kx  A)|  <  1  and  this  can  always  be  insured  by 
choosing  an  appropriate  Cz  in  (6).  In  other  words,  the  contour  Cz  is  chosen  such  that  the 
zeroes  of  the  function  1  —  z  cosfe  A)  lie  outside  it.  Substituting  this  into  (6)  and  making  use 
of  (8),  we  finally  arrive  at 

l  r/A 

u>i(mA,  je)  =  —  /  cosI(^A)[T1fe)  -  i0(j  -  1)  tanfe  A)T2(k,)]  e"‘'A  dkx,  (14) 

Z7t  J-jt/a 

where  ©(•)  is  the  Heaviside  step  function.  For  a  given  Yi  {kx)  and  Y 2fe),  integral  (14)  may  be 
computed  efficiently  by  the  application  of  the  inverse  fast  Fourier  transform  (iFFT)  algorithm 
[22].  However,  for  special  values  of  Yifc)  and  Y 2fe),  the  integral  may  be  evaluated  in 
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Figure  1.  Calculated  values  of  se)  for  T\(kx)  =  1,  T2(kx )  =  0. 


a  closed  form.  For  example,  with  uq(raA,0)  =  ^-<5® ,  w2{mA,  0)  =  0(=^  Tifc)  = 
1,  T2 (kx)  =  0)  and  m  and  5  even,  (14)  reduces  to  ([23],  3.631-17) 


w\(mA,  sc)  A 


1 

2s 


( 


(s  —  m)/ 2/ 


m  ^  s. 


(15) 


The  right-hand  side  of  (15)  gives  the  probability  of  finding  a  particle  at  m  in  5  steps,  given  that 
it  started  at  the  origin  at  5  =  0,  in  a  symmetric,  discrete-time,  ID  random  walk.  The  result  can 
be  obtained  directly  from  combinatorial  analysis  and  is  available  in  standard  texts  ([24],  p  75), 
([25],  p  16).  Figure  1  shows  a  plot  of  w\{mA,  sc) A  for  5  =  20,  30  and  40,  where  the  data  at 
discrete  m  has  been  connected  by  smooth  lines  for  the  sake  of  visual  clarity.  The  plots  clearly 
exhibit  the  diffusive  behavior  of  w\,  wherein  w\  spreads  out  in  space  with  a  diminishing  peak 
value  as  5  increases.  Using  the  identity  J2m=-oo  exp(zbirax)  =  2jr8(x),  —7t  ^  x  ^  tt,  where 
<5(-)  is  the  delta  function,  it  may  be  easily  verified  from  (14)  that  Ylm=-o o  w>i(mA,sc)A  =  1. 
Also  note  that  w\  >  0.  Hence  w\  A  behaves  like  a  true  probability  mass  function. 

We  are  also  interested  in  the  continuum  limits  A^0,e— >-0,ra— >-oo,  and  s  oo 
such  that  A2 /2c  =  D  >  0,mA  -a-  x,sc  — >►  t.  Using  the  results  lim  a^o  [cos^fe  A)]  = 

exp  (— k2Dt ),  limA^o  [cos5(/:x  A)  tan(&x  A)]  =  0  in  (14),  we  arrive  at 
1  C°° 

Wi  (x,  t)  =  —  Ti  (kx)  t~k*Dt  e'kxX  d kx .  (16) 

2^  J —oo 


This  is  the  well-known  solution  of  the  diffusion  equation  dw\/dt  =  Dd2w\/dx2  in  an 
unbounded  medium  with  an  initial  spectral  content  Yife)  (see,  for  example,  [26]).  For  an 
impulsive  initial  condition,  Yife)  =  1,  and  one  recovers  the  Green’s  function  w\(x,  t )  = 
exp(— x2 /ADt ) /V 47r Dt.  The  function  w\(mA,  sc)  given  in  equation  (14)  is  the  discrete 
version  of  w\(x,  t)  and  is  seen  to  depend  not  only  on  Yi (kx),  but  also  on  T2(kx).  The  latter 
contribution  arises  entirely  from  the  discrete  nature  of  space  and  vanishes  in  the  continuum 
limit.  To  summarize,  the  quantity  uq(raA,  56)  A  that  describes  the  probability  that  a  particle 
leaves  (m A,  sc)  in  either  direction  and  in  any  spin  state  describes  the  diffusion  process  for  a 
symmetric  4RW. 
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4.2.  Wave-like  behaviour 


The  wave-like  behavior  of  the  particle  motion  is  governed  by  the  discrete  functions  q\  and  q2. 
The  governing  equations  in  this  case  are  repeated  below  from  (4): 


Our  objective  here  is  to  derive  closed-form  expressions  for  the  transitional  probabilities  q\ 
and  q2.  Let  Qj(kx,z )  be  the  C  transforms  of  qjimA,  se),  and  let  Tj(kx)  be  the  Fourier 
transforms  of  the  initial  distribution  qj(mA,0),  j  =  1,  2.  On  applying  the  C  transform  to 
(17)  and  making  use  of  properties  (9)  and  (10)  and  carrying  out  some  algebraic  manipulations, 
we  get 


'Gi(**.z)l 

1 

“  1  z  (jkx  A 

1  V2e  ' 

z  pXkx A  “ 
V2e 

Tjfe)] 

Qi(kx,z)  J 

(1  -  V2 Z  cos(kx  A)  +  Z2) 

z  ~-ikxA 

L  V2e 

1  z  fz—ikxA 

1  V2C  J 

r2(kx)\ 

To  permit  evaluation  of  the  integral  with  respect  to  z  in  the  inverse  transform,  we  need  to 
express  Q\  and  Q2  in  a  separable  form  with  respect  to  kx  and  z.  To  this  end,  we  make  use  of 
the  identity  ([23],  8.945.2) 


1 

1  —  2  tx  + 1 2 


00 

YJUn(X)t\ 

0 


(19) 


where  £/„(•)  is  the  Chebyshev  polynomial  of  the  second  kind  of  order  n,  in  (18)  to  arrive  at 


Qi(kx,z ) 


Qi(kx,z) 


Z  —ikx At-.  ^7  x  .  ^ 


vT' xAr^\ 


( cosk,A\ 


A?  e-^iM*,)  +  (1  -  -A  e-^A  j  r2ar)J . 


V2 


(20) 

(21) 


As  with  the  diffusive  case,  the  contour  Cz  in  the  inverse  transform  is  chosen  such  that  the 
zeroes  of  the  denominator  function  (1  —  \flz  cos(kx  A)  +z2)  lie  outside  it.  Equations  (20)  and 
(21)  may  be  substituted  into  the  definition  of  the  inverse  transform  (6)  and  the  integral  with 
respect  to  z  evaluated  by  making  use  of  (8).  For  reasons  that  will  become  clear  shortly,  we 
are  interested  in  the  composite  discrete  function  fc(mA,  se)  =  q2(mA,  se)  +  iq\{mA,  se), 
which  will  be  compared  directly  with  the  solution  of  the  Schrodinger  equation.  The  expression 
for  is 

1  fnlA  \  /ms  k  A  \ 

fd(mA,se)  =  —  J  1 14  (— -^=r=J  [Fz  +  irjfe)] 

(cos k  A \ 

-J=~)  [(e_iir/4r  1  (kx)  -  e^4r2fe))cosfeA) 

+  (e-^r^*,)  +e*/4r2(kx))  sinfoA)]  J  eim^A  d kx.  (22) 

As  in  section  4.1,  the  integral  in  (22)  may  be  evaluated  efficiently  by  employing  the  iFFT 
algorithm.  In  the  special  case  of  Ti  (kx)  —  0,  T2(kx)  —  K2,  a  constant,  the  expression 
provided  in  (22)  can  be  further  simplified.  Making  a  change  of  variable  y  =  cosfeA) 
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and  using  dy/^/l  -  y2  =  —dkxA,  cos(mcos  1  y)  =  Tm(y ),  where  Tm(-)  is  the  Chebyshev 
polynomial  of  the  first  kind  of  order  ra,  we  can  show  that 


xf/dimA,  sc)  A  =  — 


=  /*' _ 1 _ { 

7T  ]_  i  ^1  -  y2  l 


v  i  A 


V2 


T’ml'y)  -  (  Af 


(ti) 


X  [rm_i (y)  +  irm+i (y)]  |  dy .  (23) 

From  the  even  and  odd  properties  of  Chebyshev  polynomials,  it  can  be  deduced  that  for 
s  —  2r  and  m  —  In  —  1  (or  vice  versa),  the  integral  in  (23)  vanishes  implying  that 
V^[( 2n  —  1)A,  2 re]  =  0  in  this  special  case. 

Other  interesting  identities  can  be  derived  starting  from  (22).  Using  the  relation 
fe(l/V2)  =  jy*  [cos(jr/4)]  =  sin(s7r/4)  +  cos(stt/4),  one  can  readily  see  that 
00 

J2  fd(mA,se)A  =  e"ilrs/4  [r2(0)  +  ir^O)] .  (24) 

m=— oo 

Hence,  unlike  w\  A,  the  quantities  q\  A  and  #2  A  can  be  of  alternating  signs  and  do  not  represent 
true  probability  mass  functions. 

Ord  [11]  has  shown  that  eight  different  continuous  functions  are  embedded  into  the 
discrete  functions  q\  and  q2.  We  will  focus  on  the  continuous  function  that  would  result 
from  choosing  x  =  2nA,n  —  0,  =bl,  ±2,  . . . ,  and  t  —  8 rc,  r  —  0,  1,  2, . . . ,  in  the  discrete 
functions  q\  and  q2.  We  show  that  x\rd  satisfies  the  Schrodinger  equation  for  m  —  2 n,  s  =  8r  in 
the  limit  as  A  ^  0,  e  — 0  ,n  —>  00  ,r—>  00  such  that  A2 /2c  =  D.  The  following  identities 
[23,  27]  involving  Chebyshev  polynomials  will  be  utilized  in  subsequent  development: 


zUs-dz)  =  us(z)  -  TS(Z) 


—  Ts(z)  —  sUs-i(z) 

dz 


(*) 


:  COS 


(1  -  z2)T''(z)  -  zT'(z)  +  s2Ts(z)  =  0 


(25) 

(26) 

(27) 


f/'-'(A)  =  '/5sin(V)- 


(1  -z2)U'\z)  -3zU'(z)+s(s+2)Us(z)  =  0,  (28) 


where  a  prime  denotes  differentiation  with  respect  to  the  argument.  For  the  purpose  of 
investigating  the  continuum  limits,  we  would  like  to  cast  (22)  in  a  form  more  suitable  for 
asymptotic  analysis.  The  last  term  in  (22)  involving  Us- 1(-)  can  be  replaced  with  dTs(-)/dkx 
on  using  the  second  relation  (26)  to  yield 
sin  kx  A 


V2 


.  /  cos  kx  A  \  —Id  /  cos  kx  A  \ 

s~l  \~sr ) =  jAdk~/s  (  41  ) 


(29) 


This  term  is  then  integrated  by  parts  and  simplified  using  the  periodicity  condition  Vj  (tt/ A)  = 
Yj(—n/ A),  j  =  1,2.  A  convenient  expression  for  the  evaluation  of  x\rd{m A,  sc)  is  then 
obtained  as 

1  ri r/A 

”'n|[ri(*,)-ir2(*x)]i/. 


fd(mA,se )  =  — 


+  (1  + 


/jr/A  t 

eimkxA\{ 

-jr/A  [ 

"(I1 


r  .m~ 

r2fe)  + 

m 

1+1- 

1+  — 

L  -s  J 

5  _ 

(cos  kxA\ 

TF" ) 

M“v2 


rife) 


cos  kx  A 


1  +  i  ,  (  cos  kx  A  \  1 

+  — tr'fe)  -  ir;fe)]r,  (— J  J  cfe, 


(30) 
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which  is  also  more  amenable  to  asymptotic  analysis  than  (22).  In  the  special  case  of 
Tifc)  =  0,  T2(kx)  =  K2  and  for  m/s  ^  0  (small  spatial  locations  and  large  times)  we 
move  on  using  (25)  that 


\//d(mA,s€) 


A  r,A  jmkx  A  I"  (  cos  kx  A  \  .  cos  kx  A 

2^ J-./a  La  V2  ;  V2 


(31) 


We  now  perform  an  asymptotic  analysis  for  small  A  in  (31)  and  show  that  i/^(2nA,  8re) 
satisfies  the  Schrodinger  equation.  To  this  end,  we  note  the  following  Taylor  series  expansions 
which  are  obtained  by  making  use  of  (26)-(28): 


cos  kx  A 


k2x  A2  (kx  A)4 

— -  +  -  +  •  •  • 

2  4! 


feA)V 

4! 


x 


+  . . . 


(32) 


(33) 


(34) 


Inserting  (32)-(34)  into  (31)  and  choosing  s  =  8r,  A2  =  2 De,  s€  =  t,mA  =  x,  s  -a-  oc, 
m  — >►  oo,  A  — >  0,  6  -a-  0,  we  arrive  at  the  desired  result: 


^ d(x,t )  =  q2(x,t)  +  iqi(x,t) 


-  i Dk\t 


In 


e-i Dk*tjkxx  ^ 


k4xD2t 2 
2! 


(35) 


Equation  (35)  is  the  spectral  representation  of  the  Green’s  function  corresponding  to  the 
Schrodinger  equation  dx/r/dt  =  iDd2\ls/dx2  with  the  impulsive  initial  condition  \lr(x,t  = 
0+)  =  K28(x).  It  has  the  exact  solution 

f(x,  t )  =  -  =  e“2/4Dg.  (36) 

V 4jtiDt 

To  reinforce  to  the  reader  that  the  plots  of  the  transitional  probabilities  (qi,  q2)  do  resemble 
the  solutions  of  the  free  Schrodinger  equation,  we  show  in  figure  2 a  comparison  of  the  real, 
9^,  and  imaginary,  S,  parts  of  the  exact  solution  (36)  of  the  Schrodinger  equation  with  the 
partial  solution  (< q\ ,  q2)  of  the  4RW.  The  numerical  solutions  shown  in  the  figure  for  i jrd  are 
on  a  discrete  spacetime  (x  =  mA,t  =  sc)  and  have  been  computed  using  (31)  with  the 
iFFT  algorithm  [22]  with  size  s  =  213  =  8192.  It  is  seen  that  the  4RW  produces  solutions 
of  oscillatory  type  with  both  positive  and  negative  excursions  for  the  expectations  q\  and  q2 , 
which  are  in  excellent  agreement  with  the  analytical  results  for  small  ™ .  This  is  in  contrast  to 
the  quantity  w\  shown  in  figure  1,  which,  behaving  like  the  solution  of  the  diffusion  equation, 
decays  exponentially  in  space  and  always  remains  positive. 
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Figure  2.  Comparison  of  the  exact  solution  of  Schrodinger  equation  with  the  discrete 
solution  of  a  4RW  for  an  impulsive  initial  condition.  ( a )  q2(mA,  se),di{\lr(mA,  se)}  and  ( b ) 
q\  (mA,  se),  se)}. 


5.  Summary 

By  considering  a  multistate  random  walk  on  a  discrete  lattice,  expressions  have  been  derived 
for  the  various  transitional  probabilities  using  the  concept  of  generating  functions.  A  2D 
transform  involving  Fourier  transformation  in  space  and  the  z -transformation  in  time  is 
employed  to  accomplish  this.  The  transitional  probabilities  governing  particle  motion  are 
expressed  in  terms  of  integrals  involving  trigonometric  functions  in  the  case  of  the  diffusion 
equation,  and  involving  Chebyshev  polynomials  of  the  first  and  second  kinds  in  the  case 
of  the  Schrodinger  equation.  Closed-form  expressions  have  been  given  for  particular  cases 
of  the  initial  conditions.  The  continuum  limits  of  the  diffusion  equation  and  Schrodinger 
equation  have  been  shown  to  follow  directly  from  these  transitional  probabilities  through 
the  performance  of  appropriate  asymptotic  analysis.  The  present  analysis  consisting  of  the 
z -transform  along  the  time  axis  and  Fourier  transform  along  the  spatial  axis  is  most  suitable 
for  studying  linear  difference  equations  with  constant  coefficients.  In  the  4RW  model,  this 
would  correspond  to  the  free  Schrodinger  equation.  The  important  extension  of  this  analysis  to 
higher  dimensions  is  worth  exploring  and  would  be  taken  up  in  the  future.  The  incorporation 
of  a  smooth  potential  field  in  the  Schrodinger  equation  into  the  4RW  model  has  already  been 
addressed  by  Ord  in  [19]  and  the  study  of  its  transitional  probabilities  will  be  taken  up  in  a 
separate  paper  using  a  different  approach. 
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Abstract — Transparent  boundary  condition  in  a  2D-space  is  pre¬ 
sented  for  the  four-state  random  walk  (4RW)  model  that  is  used 
in  treating  the  standard  parabolic  equation  by  stochastic  methods. 
The  boundary  condition  is  exact  for  the  discrete  4RW  model,  is 
of  explicit  type,  and  relates  the  field  in  the  spectral  domain  at  the 
boundary  point  in  terms  of  the  field  at  a  previous  interior  point 
via  a  spectral  transfer  function.  In  the  spatial  domain,  the  domain 
of  influence  for  the  boundary  condition  is  directly  proportional  to 
the  “time”  elapsed.  By  performing  various  approximations  to  the 
transfer  function,  several  approximate  absorbing  boundary  condi¬ 
tions  can  be  derived  that  have  much  more  limited  domain  of  influ¬ 
ence. 

Index  Terms — Generating  function,  parabolic  equation,  random 
walk,  Schrodinger  equation,  transform  methods,  transparent 
boundary  condition. 


I.  Introduction 


PARABOLIC  EQUATION  (PE)  is  used  widely  in  several 
areas — including  radiowave  propagation,  underwater 
acoustics,  fiber  optic  propagation,  etc. — to  study  propagation  of 
unidirectional  waves  in  a  variety  of  environments  [5],  [8],  [11]. 
Its  form  is  very  similar  to  the  Schrodinger  equation  in  which 
the  wavefunction,  temporal  variable,  and  the  potential  take  the 
respective  places  of  the  reduced  field,  the  axial  coordinate,  and 
the  modified  refractive  index  of  the  PE.  When  solving  the  PE  or 
Schrodinger  equation  in  unbounded  regions  by  differential  or 
difference  equation  methods,  one  has  to  first  restrict  the  region 
to  a  finite  size  and  simulate  an  unbounded  region  by  imposing 
boundary  conditions  on  the  truncated  domain.  If  these  boundary 
conditions  exactly  simulate  an  unbounded  exterior  region  as  far 
as  the  interior  problem  is  concerned,  they  are  termed  as  exact  or 
transparent  boundary  conditions.  Transparent  boundary  condi¬ 
tions  can  be  devised  for  the  original  partial  differential  equation 
or  to  its  discretized  version,  and  this  has  been  done  by  several 
works  in  the  past:  [2],  [4],  [9],  [10],  [13],  and  [1].  ([1]  and  [4] 
contain  many  works  outside  the  traditional  electromagnetic  area 
where  such  development  took  place.)  For  a  number  of  reasons, 
there  is  an  interest  to  treat  the  PE  by  stochastic  techniques, 
wherein  the  desired  field  is  obtained  by  performing  averages 
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over  several  realizations  of  Brownian  motion  or  its  discrete 
counterpart  of  random  walk  traversing  the  problem  domain 
[6].  A  direct  application  of  stochastic  techniques  to  the  PE  as 
done  in  [3]  necessitates  analytical  continuation  of  fields  and 
boundary  data,  only  possible  for  very  limited  problems,  that  is 
highly  undesirable.  The  four-state  random  walk  (4RW)  model, 
originally  developed  by  Ord  [12]  to  provide  a  macroscopic 
model  for  the  Schrodinger  equation,  alleviates  this  limitation 
and  results  in  the  same  second-order  accurate  discretization 
scheme  for  the  PE  as  the  Crank-Nicolson  scheme.  Numerical 
schemes  are  currently  being  developed  using  the  4RW  model 
to  study  wireless  propagation  problems  in  open  domains,  and 
it  is  the  purpose  of  this  letter  to  derive  transparent  boundary 
conditions  for  this  model  by  considering  2D  propagation  in  a  1 
+  1 -dimensional  space.  Detailed  numerical  solutions  for  prop¬ 
agation  over  terrain  and  rough  surfaces,  and  in  inhomogeneous 
atmosphere  by  the  4RW  model  will  be  presented  in  a  future 
paper.  In  the  next  section,  we  begin  with  a  brief  introduction 
of  the  4RW  model,  followed  by  the  development  of  the  exact 
boundary  condition. 


II.  Transparent  Boundary  Condition  for  4RW 


Consider  the  standard  parabolic  equation  in  the  reduced  field 
variable  ^(x,  t): 


d ^ 
dt 


i  d 2/0 
2&o  dx2 


+  ikoVip 


(1) 


which  constitutes  a  narrow  angle  approximation  to  the 
Helmholtz  equation  for  time-harmonic  waves  (e~'lk°T  time 
convention  in  the  normalized  time  variable  r)  with  wavenumber 
ko  =  2tt/X  propagating  in  the  t-x  plane,  where  A  is  the  op¬ 
erating  wavelength  and  V(x,t)  is  a  complex- valued  function 
proportional  to  the  modified  refractive  index  (or  potential 
for  the  Schrodinger  equation)  of  the  medium.  Because  this 
work  relies  heavily  on  the  random  walk  literature  pertaining 
to  the  Schrodinger  equation,  we  shall  refer  to  t  as  the  “time” 
coordinate  even  though  it  represents  the  axial  space  coordi¬ 
nate  in  wave  propagation  problems.  The  variable  x  represents 
the  lateral  spatial  coordinate  and  i  =  \/“T-  Presence  of  the 
first-order  derivative  in  time  enables  the  solution  of  (1)  to 
proceed  with  the  specification  of  initial  conditions  at  t  =  0.  It 
is  assumed  that  all  sources  at  initial  time  are  confined  to  the 
region  x  <  0  (Region  I)  and  that  the  potential  V  assumes  a 
constant  value  V(x,t)  =  Vc  for  x  >  0  (Region  II).  Our  interest 
is  to  derive  an  artificial  boundary  condition  imposed  at  x  =  0 
for  a  discretized  version  of  (1)  that  simulates  free  space  as  far  as 
Region  I  is  concerned.  For  simplicity  and  no  loss  in  generality, 
we  first  consider  the  case  of  Vc  =  0  and  discuss  the  solution 
for  the  nonzero  case  at  the  end.  For  treatment  by  stochastic 
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methods,  (1)  is  discretized  as  per  the  4RW  model,  which, 
like  the  Crank-Nicolson  finite-difference  scheme,  results  in  a 
second-order  accurate  (in  space  and  time)  scheme  [12].  For  a 
particle  moving  on  a  discrete  lattice  and  subject  to  random  col¬ 
lisions,  the  transitional  probabilities  in  the  4RW  at  the  discrete 
space-time  point  {x  =  mA,  t  =  se)  are  of  the  form  [7],  [12] 


(2) 


where  qi(mA,se)  =  2S//2[pi(raA,  se)  —  ps{mA^se)], 

q2(mA,se)  =  2s/2[p2(mA,  se)  —  p^mA^se)],  and 

pM(mA,  se)  A,  ft  =  1, . . . ,  4,  is  the  probability  that  a  particle  is 
in  state  /i  at  the  space-time  point  (mA,  se),  ra  =  0,  ±1,  ±2, . . ., 
s  =  0, 1, . . ..  The  particle  changes  its  direction  of  motion  with 
every  collision  but  changes  its  parity  or  spin  at  every  other 
collision.  The  combination  of  two  directions  of  motion  and  two 
states  of  spin  constitutes  the  four  states  in  the  model.  It  has 
been  shown  in  [12]  that  such  a  4RW  encompasses  the  diffusion 
as  well  as  Schrodinger  equations  and  that  =  g2  +  i(]i  is 
the  Schrodinger  wave  function  in  the  discrete  case.  The  rela¬ 
tionship  between  the  temporal  and  spatial  steps  of  the  form 
6  =  koA 2  implied  in  the  4RW  model  assures  the  aforemen¬ 
tioned  second-order  accuracy  in  space.  For  waves  propagating 
within  ±15°  about  the  axial  direction,  the  lateral  space  step  is 
usually  restricted  to  A  <  2 A  [8],  resulting  in  e  <  87rA,  values 
which  are  consistent  with  those  used  in  the  normal  finite-differ¬ 
ence  or  Fourier  split-step  schemes  [8]  used  to  solve  the  PE. 

The  operators  Ex  and  Et  in  (2)  are,  respectively,  the 
spatial  and  temporal  advancing  operators  and  are  defined 
mathematically  as  E^p^mA,  se)  =  pfl  [{m  =b  1)A,  se]  and 
EtPfi{mA^se)  =  p^[mA,  (s  +  l)e].  It  is  assumed  in  (2)  that 
the  random  walk  is  nonpersistent  and  the  probabilities  that  a 
particle  maintains  or  changes  its  direction  at  the  next  time  step 
are  the  same.  For  brevity  we  will  denote  ?/;(mA,  se)  as  simply 
s).  As  in  our  previous  studies,  we  will  assume  that  the 
number  of  right-going  particles  is  the  same  as  those  going  to  the 
left  at  time  t  =  0  and  focus  on  deriving  transparent  boundary 
conditions  for  the  wave  function  $(m,s).  The  boundary  at 
m  —  0  separates  Region  I  from  Region  II  and  it  is  further 
assumed  for  convenience  that  0)  =  0  for  m  >  —2. 
Consider  a  discrete  function  v(m,  s)  and  its  temporal  transform 
v(rri.  z)  as  defined  in  [7] 


v{m,  z)  =  v(m,  s)zs  =  Zv.  (3) 

,s=0 

The  quantity  v(m,z)  may  be  thought  of  as  the  discrete  version 
of  the  Laplace  transform  of  v.  The  inverse  relation  is  defined  as 

s)  =  /  V~^rdz  =  z~li}  (4) 

6* 

where  Cz  is  a  closed  contour  around  the  origin  in  the  com¬ 
plex  z -plane  that  encloses  only  the  singularities  at  the  origin. 
The  symbols  Z  and  Z~x  denote  the  temporal  transform  and  its 
inverse  respectively.  The  transformed  variable  t)(m,  z)  is  also 
referred  to  as  the  generating  function  within  the  random  walk 
community  [14].  If  v\  (rn.  s)  and  vs  (rn,  s)  are  two  temporal  (and 


causal)  functions  with  respective  generating  functions  v\  (m,  z) 
and  V2 (m,  z),  the  following  shift  and  convolution  properties,  re¬ 
spectively,  can  be  readily  established: 

2[Etv i]  =  z_1  [Di  (m,z)  -^i(to,0)]  (5) 

S 

z-'lvihl  =  E  v\(m,  s  —  k)v2(m,  k).  (6) 

k= o 

Applying  the  temporal  transform  to  (2)  and  making  use  of  the 
shift  property  (5),  and  carrying  out  some  algebraic  manipula¬ 
tions,  we  arrive  at  the  following  equations  for  the  generating 
function  in  terms  of  initial  data: 

[1  —  V2zDx  +  z2\f(m,  z)  —  *$&(ra,  0) 

[(E*l+iEx) ?i(m>°)]  (?) 

where  Dx  =  (Ex  +  Ef1)/ 2  is  the  discrete  averaging  operator. 
The  right-hand  side  in  (7)  is  identically  zero  for  Region  II,  and 
a  knowledge  of  the  general  solution  in  this  region  permits  the 
derivation  of  transparent  boundary  condition.  A  general  solution 
of  (7)  in  Region  II  will  be  of  the  form  'f(rriy  z)  =  A[0(z)]rn  with 
A  constant.  The  function  v(z)  is  then  determined  from  (7)  as 

Hz)  =  -T-[i  +  ^2  +  yi  +  z4].  (8) 

V2  z 

We  are  interested  in  the  solution  that  remains  bounded  for  large 
m  in  Region  II,  which  imposes  the  condition  |  v(z)  |  <  1.  Due  to 
the  presence  of  the  radical  sign,  the  function  defined  by  (8)  will 
be  multiple-valued,  resulting  in  two  branches  of  the  solution. 
These  two  solutions  will  cross  over  in  the  complex  plane  unless 
care  is  taken  to  define  each  analytically  by  considering  a  double- 
sheeted  complex  plane.  The  two  sheets  intersect  at  the  branch 
cuts  defined  by  the  curve  |  v{z)  \  —  1  so  as  to  separate  the  proper 
solution  \v(z)\  <  1  from  the  improper  one  \v{z)\  >  1.  It  is 
more  illuminating  to  first  study  the  branch  cuts  in  a  transformed 
complex  plane  (  =  (r  +  iQ  =  z2,  in  which  the  branch  points 
occur  at  (  =  ±i.  In  the  z-plane,  the  branch  points  occur  at  Zk  = 
exp[7(2&  —  l)7r/4],  k  =  1, . . . ,  4.  With  these  branch  cuts,  the 
top  sheet  of  the  complex  z-plane  is  chosen  such  that  y/l  +  z4  = 

\/nt=i(z  —  zk)  =  —  1  at  2  =  0.  A  straightforward  but  tedious 
analysis  reveals  that  the  branch  cut  in  the  £-plane  is  given  by  the 
semicircle  \(\  =  1,  with  the  real  part  £r  <  0.  The  branch  cuts  in 
the  (-  and  jz-  planes  are  shown  in  Fig.  1.  We  denote  the  proper 
branch  by  the  symbol  v\  (z)  and  the  other  branch  by  the  symbol 
i>2  (z).  For  the  purpose  of  numerical  calculations,  one  may  set 
z)i  (z)  =  [1  +  z2  =p  x/l  +  z4]/\/2z.  Clearly  ki(z)i>2(z)  =  1. 

The  proper  solution  in  Region  II  is  now  given  by  fj(m,z)  = 
v4[z>i(z)]m,  which  implies  that 

^(o  ,z)  =  z>i(4^(-l,4.  (9) 

Equation  (9)  is  the  frequency  domain  representation  of  the  trans¬ 
parent  boundary  condition  that  simulates  free  space  for  the  in¬ 
terior  region  m  <  0.  The  function  k{z)  may  be  interpreted  as 
a  transfer  function,  which  relates  the  input  and  output  functions 
/0(— 1  ,z)  and  z),  respectively.  For  the  purpose  of  numer¬ 
ical  implementations,  it  is  necessary  to  prescribe  the  boundary 
condition  in  the  temporal  domain.  From  Fig.  1,  it  is  clear  that 
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the  function  v\  (z)  is  analytic  near  the  origin.  The  inverse  trans¬ 
form  of  v\(z)  involves  integration  around  the  closed  contour 
Cz  and  is  facilitated  by  performing  a  Taylor  series  expansion 
of  Vl  +  near  the  origin.  The  infinite  series  representation  of 
z> i  (z)  can  then  be  easily  determined  as 


+  ^  arzAr+7 

r= 0 


(10) 


where 


ar 


(— l)r(2r  +  1)! 

4r+ir\(r  +  2)l\/2 


0(r  3/2),  r  — ►  oo. 


Making  use  of  the  relation 


(ID 


Fig.  1.  Branch  cuts  in  the  complex  (  =  z2  and  2  planes  dehning  |  (  z )  |  =  1. 

The  branch  points  in  the  (-plane  occur  at  (i  =  i  and  (2  —  —i  and  those  in  the 
z -plane  at  zk  —  eA2^-1)^/4^  fc  —  17 .  .  .  ?  4.  Branch  cuts  are  shown  as  wiggly 
lines. 


where  6k  is  the  Kronecker’s  delta,  we  arrive  at  the  corre¬ 
sponding  temporal  function 


Vl  (s)  =  4=6* - 1 

W  V2  2V2 


83  +  ar8f+7. 

r= 0 


(12) 


Finally,  a  “time”-domain  transparent  boundary  condition  is  ob¬ 
tained  by  using  the  convolution  property  (6)  in  (9) 


m  s)  =  1,  a  -  1)  -  2  S  -  3) 


R 


+  ^  ar^{— 1,  s  —  7  —  4 r)  (13) 


r=0 


where  —  7)/4j  is  the  interger- valued  floor  function  of 

(. s  —  7) /4.  Obviously,  the  summation  term  will  only  exist  when 
s  >  7.  Note  that  the  presence  of  a  term  zk  in  the  series  rep¬ 
resentation  of  v\{z)  results  in  a  delayed  field  ^(  — 1,  s  —  k)  in 
the  temporal  representation  of  the  boundary  condition.  Equation 
(13),  consisting  of  a  finite  series,  is  the  exact  discrete  boundary 
condition  for  terminating  the  computational  domain  of  the  4RW 
model  and  is  based  on  an  infinite  series  representation  of  the 
transfer  function  in  a  region  analytic  about  the  origin.  Further¬ 
more,  it  is  of  explicit  type  wherein  the  field  at  (0,  s)  is  expressed 
in  terms  of  the  historical  field  at  (  —  1,  s  —  k),  k  >  0.  It  is  also 
seen  that  the  number  of  terms  in  (13)  will  increase  linearly  with 
s.  To  reinforce  to  the  reader  its  accuracy,  we  place  initial  sources 
</i(ra, 0)A  =  Ti^°,  q2(m,0)A  =  T2<5™°  with  constant  Ti 
and  T 2  and  compare  the  exact  solution  of  the  field  (0,  s)  as  de¬ 
termined  from  the  expressions  in  [7]  with  the  solution  ^(0,  5) 
obtained  by  substituting  the  exact  solution  into  the  right  hand 
side  of  (13).  Note  that  the  earliest  “time”  at  which  an  impulse 
placed  initially  at  mo  =  —  28  will  arrive  at  m  =  0  is  s  =  28. 
The  two  solutions  plotted  in  Fig.  2  for  s  >  32  are  seen  to  be 
virtually  on  top  of  each  other. 

It  is  possible  to  derive  a  hierarchy  of  approximate  boundary 
conditions  by  considering  various  approximations  to  the 
transfer  function  For  example,  one  may  perform  a 

rational  approximation  to  the  function  and  obtain  suc¬ 

cessively  higher  order  boundary  conditions.  A  zeroth  order 
boundary  condition  is  obtained  by  using  v\(z)  —  zj\f2  and 


Fig.  2.  True  field  in  unbounded  medium  and  that  obtained  from  (13).  (a)  Real 
part  of  s).  (b)  Imaginary  part  of  r(0.  .s  ). 


reads  5)  =  /0(— 1,  s  —  1).  A  better  approximation  uses 

v\{z)  —  y/2z/{2  +  z2)  and  results  in  the  boundary  condition 
2-0(0,  s)  +  0(0,3  —  2)  =  a/20(— l,s  —  1).  Approximations 
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in  terms  of  special  functions  such  as  Legendre  polynomials  are 
also  possible. 

It  is  straightforward  to  handle  the  case  with  nonzero 
Vc  by  writing  the  solution  <j){x,t)  of  (1)  as  = 

/)  exp  (ikoVcL),  resulting  in  s)  =  s) 

with  r*i  =  ex.p(ikoeVc) .  The  boundary  condition  (9)  con¬ 
tinues  to  be  valid  for  the  quantity  The  spectral  domain 
boundary  condition  in  terms  of  the  actual  field  becomes 
</>( 0,  z)  =  1,  z).  In  the  temporal  domain  we  get  the 

corresponding  equation 


V(0, S)  =  ^(-1, *  - 1)  -  5  -  3) 


+  ^2  arr\r+1 </>(-!,  s-7  -  4 r).  (14) 


Obviously,  (14)  reduces  to  (13)  for  Vc  =  0  (or  n  =  1). 


III.  Conclusion 

Transparent  boundary  condition  has  been  presented  for  the 
4RW  model  used  to  represent  a  PE  in  stochastic  methods.  The 
boundary  condition,  given  in  the  form  of  a  finite  series,  involves 
specification  of  the  field  at  a  boundary  point  in  terms  of  the  his¬ 
torical  field  at  immediate  interior  points.  Several  approximate 
boundary  conditions  can  be  derived  from  it  based  on  a  rational 
or  other  approximations  to  the  transfer  function. 
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An  ADI-PE  Approach  for  Modeling  Radio 
Transmission  Loss  in  Tunnels 
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Abstract — Alternate  direction  implicit  (ADI)  method  is  used 
to  study  radio  wave  propagation  in  tunnels  using  the  parabolic 
equation  (PE).  We  formulate  the  ADI  technique  for  use  in  tunnels 
with  rectangular,  circular  and  arched  cross  sections  and  with  lossy 
walls.  The  electrical  parameters  of  the  lossy  walls  are  charac¬ 
terized  by  an  equivalent  surface  impedance.  A  vector  PE  is  also 
formulated  for  use  in  tunnels  with  lossy  walls.  It  is  shown  that  the 
ADI  is  more  computationally  efficient  than  the  Crank  Nicolson 
method.  However,  boundary  conditions  become  more  difficult  to 
model.  The  boundary  conditions  of  the  ADI  intermediate  planes 
are  given  the  same  boundary  conditions  as  the  physical  plane  and 
the  overall  accuracy  is  reduced.  Also,  when  implementing  the  ADI 
in  tunnels  with  circular  cross  sections  the  order  at  which  the  line 
by  line  decomposition  occurs  becomes  important.  To  validate  the 
ADI-PE,  we  show  simulation  results  for  tunnel  test  cases  with 
known  analytical  solutions.  Furthermore,  the  ADI-PE  is  used  to 
simulate  real  tunnels  in  order  to  compare  with  experimental  data. 
It  is  shown  that  the  PE  models  the  electric  fields  most  accurately 
in  real  tunnels  at  large  distances,  where  the  lower  order  modes 
dominate. 

Index  Terms — Alternate  direction  implicit  (ADI),  parabolic 
equation,  radio  wave  propagation  in  tunnels. 


I.  Introduction 

THE  growth  of  mobile  communication  systems  in  recent 
years  has  led  to  increasing  interest  in  the  research  of  radio 
wave  propagation  in  tunnels.  The  methods  that  have  been  tra¬ 
ditionally  used  to  model  radiowave  propagation  in  tunnels  are 
modal  analysis,  geometrical  optics  and  the  parabolic  equation 
(PE)  approximation. 

The  modal  analysis  method  and  the  geometrical  optics 
method  both  have  unacceptable  limitations  and  cannot  always 
be  used  to  solve  for  the  fields  in  real  tunnels  [3]— [5] .  The  modal 
approach  uses  a  simple  field  representation  comprised  of  one 
dominate  mode  and  an  infinite  number  of  higher  order  modes. 
This  representation  is  not  always  applicable  to  real  tunnel 
environments  with  arbitrary  geometries.  Also,  it  is  difficult  to 
determine  eigenfunctions  for  real  tunnels  [4].  The  geometrical 
optics  method  approximates  the  field  as  rays  reflected  along 
the  direct  path  of  propagation.  The  geometrical  optics  method 
becomes  difficult  to  use  for  problems  at  long  range  because 
of  the  large  number  of  reflected  waves,  and  the  method  fails 
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Fig.  1.  Tunnels  with  (a)  rectangular,  (b)  circular  and  (c)  arched  profiles. 

completely  in  caustic  regions.  By  comparison,  the  standard 
parabolic  equation  has  been  shown  to  accurately  model  electro¬ 
magnetic  fields  in  tunnels  for  waves  which  travel  within  ±15° 
to  the  axis  of  propagation  [2] . 

Research  has  been  done  on  the  parabolic  approximation  of 
radiowave  propagation  in  tunnels  by  Popov  and  Zhu  [4]  and 
Noori  et  al.  [5].  The  traditional  numerical  technique  used  in  pre¬ 
vious  works  is  the  Crank  Nicolson  method.  The  popularity  of 
the  Crank  Nicolson  method  is  due  to  the  fact  that  it  is  stable  for 
any  discretization  in  the  transverse  plane  or  along  the  propaga¬ 
tion  axis  [2].  The  major  limitation  of  the  CN  method  is  that  it 
requires  the  solution  of  sets  of  simultaneous  equations  that  may 
become  too  large  to  efficiently  solve  for  problems  with  dense 
meshes.  The  alternate  direction  implicit  (ADI)  technique  was 
developed  to  address  the  problem  of  computational  efficiency. 
The  ADI  method  has  been  used  extensively  in  previous  works  to 
solve  the  PE  in  the  form  of  the  heat  transfer  equation  [19]— [22] 
and  in  electromagnetics  to  directly  solve  Maxwell’s  equation 
using  the  FDTD  method  [23]— [25]. 

In  this  paper  we  use  PE  together  with  the  ADI  technique  to 
study  transmission  loss  inside  tunnels.  In  Section  II,  we  show 
how  the  PE  can  be  used  to  solve  for  fields  in  tunnels  with  per¬ 
fectly  conducting  walls.  In  this  case,  we  can  use  a  scalar  PE  and 
enforce  simple  Dirichlet  or  Neumann  boundary  conditions  on 
the  wall.  To  model  waves  in  tunnels  with  lossy  walls,  we  must 
use  a  vector  parabolic  equation,  as  discussed  in  Section  III.  A 
brief  introduction  of  the  CN  method  is  given  in  Section  IV.  As 
discussed  in  Section  V,  the  ADI  maintains  the  unconditional  sta¬ 
bility  of  the  CN  method  while  reducing  compuational  labor.  We 
develop  the  ADI  to  deal  with  straight  tunnels  with  rectangular, 
circular  and  arched  profiles,  as  shown  in  Fig.  1. 

For  validation  purposes,  in  Section  VI,  we  examine  the  ac¬ 
curacy  of  the  scalar  parabolic  equation  when  used  to  solve  for 
fields  in  (a)  square  (b)  circular  (c)  semi-circular  cylindrical  PEC 
tunnels.  Then  we  examine  the  accuracy  of  the  vector  parabolic 
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equation  when  used  to  solve  for  fields  in  (d)  square  and  (e)  cir¬ 
cular  cylindrical  tunnels  with  impedance  walls  in  Section  VII. 
Finally,  in  Section  VIII  we  compare  the  ADI  simulation  results 
with  published  measured  data  for  the  case  of  (f)  arched  cylin¬ 
drical  tunnels. 

II.  Scalar  Parabolic  Equation 

The  parabolic  equation  is  an  approximation  of  the  Helmholtz 
equation  that  assumes  wave  propagation  in  predominantly  one 
direction  [2].  The  PE  can  be  obtained  from  the  Helmholtz 
equation 

/  qz  qz  \ 

\^  +  W2  +  d^  +  k2°)^{X'mZ)  =  {)  (1) 

where  is  either  the  electric  or  magnetic  scalar  potential  and 
kQ  is  the  free  space  wave  number.  We  choose  a  solution  for  the 
scalar  potential  with  a  fast  varying  phase  term  in  the  2  direction 
factored  out 


^(rc,  y,  z)  =  u(x,  y,  z)e  jk°z.  (2) 

where  u(x ,  y ,  2)  is  the  reduced  plane  wave  solution.  The  re¬ 
duced  plane  wave  is  slowly  varying  along  the  2  direction  and 
we  can  make  the  following  assumption, 


d2u 
d  z2 


<kQ 


du 

dz 


(3) 


Physically,  the  reduced  plane  wave  represents  the  lower  order 
modes  which  travel  predominantly  along  the  2-axis.  Substi¬ 
tuting  the  plane  wave  solution  into  (1),  the  standard  PE  follows: 


du  1  /  d2  d2  \ 

dz  2 jkQ  \  Ox2  +  dy2  )  U' 

The  PE  reduces  the  2nd  order  2  derivative  into  a  first  order 
derivative  and,  as  a  result,  the  PE  is  easier  to  implement  numer¬ 
ically  using  finite  differences.  As  shown  in  [7],  the  relative  error 
epression  of  the  PE  can  be  shown  to  be 


^3  DPE  =  1  —  cos(0)  — 


(5) 


where  the  angle  0  is  the  angle  of  propagation  measured  from 
the  2  axis.  The  level  of  accuracy  will  determine  the  range  of 
angles  over  which  the  PE  is  valid.  For  this  paper,  we  adopt  the 
traditional  angle  range  of  ±15° .  For  the  angle  of  15°  the  relative 
error  is  0.00058.  The  PE  is  only  valid  for  waves  traveling  in  the 
positive  2  direction  at  angles  within  ±15°  and  back  scattered 
waves  are  ignored  [2] . 

In  Section  VI,  we  use  the  PE  to  solve  for  fields  in  PEC  waveg- 
iudes.  For  these  examples,  we  choose  our  initial  fields  so  that 
only  modes  propagating  within  ±15°  are  significant.  The  scalar 
PE  can  only  be  used  for  tunnels  with  PEC  walls.  In  Section  III, 
we  discuss  the  modifications  that  must  be  made  to  the  PE  to  deal 
with  tunnels  with  non-PEC  walls. 


III.  Vector  Parabolic  Equation 

The  scalar  PE  is  only  valid  in  cases  where  the  transverse  mag¬ 
netic  and  electric  fields  can  propagate  independantly  of  each 
other.  This  is  the  case  when  studying  perfect  electrically  con¬ 
ducting  (PEC)  waveguides.  If  Dirichlet  boundary  conditions  are 
enforced  on  the  waveguide  wall,  the  propagating  field  will  be  of 
the  TM  type,  and  if  Neumann  boundary  conditions  are  enforced, 
then  the  field  will  be  of  the  TE  type. 

In  realistic  tunnels,  the  walls  are  typically  made  up  of  some 
lossy  material  such  as  concrete  and  rock.  The  PE  approxima¬ 
tion  has  been  shown  to  model  realistic  tunnels  very  accurately 
over  large  distances  [4] .  This  is  because  higher  order  modes  are 
more  severely  attenuated  over  large  distances  due  to  multiple  re¬ 
flections.  Lower  order  modes,  which  propagate  at  small  angles 
relative  to  the  axis,  make  the  most  significant  contribution  to  the 
field  over  large  distances. 

Also,  in  the  case  of  realistic  tunnels,  there  will  exist,  in  gen¬ 
eral,  coupling  between  the  TM  and  TE  modes.  Unlike  the  PEC 
case,  we  cannot  solve  the  TE  and  TM  modes  independently.  We 
must  now  develop  a  vector  PE  that  will  take  coupling  effects 
into  account.  Formulation  of  a  vector  PE  for  tunnels  with  lossy 
walls  and  a  radius  of  curvature,  p,  was  done  by  Popov  and  Zhu 
[4].  For  this  paper,  we  consider  only  straight  tunnels  with  infi¬ 
nite  radius  of  curvature.  Although  we  will  not  repeat  the  details 
here,  we  will  briefly  discuss  the  vector  PE  derivation.  In  [4],  the 
electrical  parameters  on  the  tunnel  walls  are  approximated  by 
the  Leontovich  impedance  boundary  condition  [27] 


n  x  E  =  —  Zn  x  [n  x  77)  (6) 

where  n  is  the  outward  normal  on  the  tunnel  wall  and  Z  is  the 
relative  surface  impedance.  For  a  wall  with  relative  permittivity 
er  and  conductivity  aQ  (in  S/m),  Z  is  approximately  [18] 


6rc 


where  erc  =  er—jar,ar  =  (j0/uje0 ,  is  the  complex  permittivity 
and  relative  conductivity,  respectively.  Here  we  use  the  grazing 
angle  approximation  for  the  surface  impedance  because  the  rays 
that  make  the  significant  contributions  to  the  field  at  large  dis¬ 
tances  propagate  at  small  angles  relative  to  the  axis  of  propaga¬ 
tion. 

Starting  with  normalized  Maxwell’s  equations 

V  x  rj0h  =  V  x  H  =  jk0E 
and 


Vx7  =  —jk0r]0h  =  —jk0H  (8) 

where  y0  is  the  intrinsic  impedance,  E  is  in  V/m ,  h  is  in  A/m 
and  H  is  the  normalized  magnetic  field  and  is  in  AEt/m.  We 
can  define  a  six  component  vector 


n-  [Ex,Ey,Ez,Hx,Hy,Hz\  . 


(9) 


MARTELLY  AND  JANASWAMY:  AN  ADI-PE  APPROACH  FOR  MODELING  RADIO  TRANSMISSION  LOSS 


1761 


We  can  then  use  the  following  asymptotic  ansatz 


II  —  U  exp 


(10) 


where 


U2  n3  /77 

C  =  \  — - and  v  =  \  —  <  1  (11) 

V  p  V 9 

and  D  is  the  tunnel  width  [26] .  The  eikonal  is  a  general  second 
order  polynomial  of  v  and  the  vector  amplitude  U  is  an  asymp¬ 
totic  series  in  powers  of  v.  Directly  substituting  the  ansatz  into 
Maxwell’s  equation  and  equating  terms  of  the  same  order  of 
v ,  we  can  define  the  eikonal  and  find  the  vector  PE  [4]  for  the 
straight  waveguide 


^  dW 
2  = 
oz 


ifiW  (J2W 
dx2  +  dy2 


(12) 


where  W  is  the  attenuation  function  which  describes  the 
complex  wave  amplitude  of  the  approximate  plane  wave 
exp  {—jkGz)  propagating  along  the  waveguide  axis.  W  is 
related  to  the  electric  field  transverse  to  the  direction  of  prop¬ 
agation  by 


£tra„s  =  (Ex,  Eyf  =  We-*h*  (13) 

Substituting  the  ansatz  into  the  IBC  [see  (6)],  and  equating  like 
terms,  we  get  [4],  [28] 


j  =  =  =  dW 

| boundary  —  ^  T0G0T0  (boundary 


(14) 


where  T0  and  G0  are  matrices  defined  by 


T  — 

-L  n  — 


nT 


-nx 


and  G0  = 


I  0 
0  z 


(15) 


where  nx  and  ny  are  the  unit  normal  vectors.  The  impedance 
boundary  condition  describes  the  effects  of  grazing  angle  reflec¬ 
tion,  selective  mode  absorption  and  depolarization  on  the  wave¬ 
guide  walls  [4] .  The  Crank-Nicolson  method  has  been  widely 
used  to  model  both  the  scalar  and  vector  PE.  In  Section  IV,  we 
will  briefly  introduce  the  Crank-Nicolson  method  as  it  will  di¬ 
rectly  lead  to  the  ADI  technique. 


IV.  Crank-Nicolson  Method 


The  Crank-Nicolson  method  has  been  widely  used  to  model 
electromagnetic  wave  propagation  in  tunnels  [3]-[5].  The  main 
feature  of  the  Crank-Nicolson  method  is  its  stability  for  any 
mesh  ratio  [11],  [17].  The  enumeration  scheme  of  the  marching 
planes  for  the  Crank-Nicolson  method  is  shown  in  Fig.  2.  The 
coordinates  of  the  mesh  points  on  the  x,  y  and  z  axes  are  denoted 
by  mAx ,  l Ay  and  nAz,  respectively.  The  field  points  on  the  n 
and  the  n  + 1  marching  planes  are  explicitly  solved,  but  the  field 
points  on  the  intermediate  plane  are  never  explicitly  calculated. 
The  Crank-Nicolson  descretization  of  the  standard  PE  is  given 
by  (16),  shown  at  the  bottom  of  the  page,  which  can  be  re-written 
as 


4 jkQ 


Jy-s 

4jk0  v 


n+1 


=  1  + 


4 jk( 


-fix  + 


—^-6 
4jk0  y 


(17) 


where  rx  =  Az/Ax 2  and  ry  =  Az/Ay 2  are  the  mesh  ratios, 
and 


fix/Urn.l  — 1,1  “b  Um—l.l  (18) 

fiyUm.J  =  2  Um^l  +  Umj  —  !•  (19) 

The  Crank-Nicolson  method  is  implicit,  unconditionally  stable, 
and  second  order  accurate  in  x,  y  and  $  [9] .  The  Crank-Nicolson 
difference  scheme  is  centered  about  the  midpoint  between  the 
z  =  nAz  and  2  =  (n  +  l)Az  marching  planes.  The  field  points 
of  each  successive  marching  plane  must  be  solved  in  consecu¬ 
tive  order  at  propagation  steps  of  Az  until  the  field  at  the  desired 
range  is  solved.  The  major  limitation  of  the  Crank-Nicolson 
method  is  that  it  becomes  computationally  intensive  in  cases  that 
require  a  very  large  mesh.  As  Equation  (16)  shows,  the  size  of 
the  matrix  generated  for  a  problem  with  Nx  mesh  points  along 
the  x  axis  and  Ny ,  mesh  points  along  the  y  axis  will  be  — 4 

(in  Fig.  2,  Nx  =  Ny  =  4).  The  elements  of  the  matrix  will 
not  be  in  tridiagonal  form.  Furthermore,  even  though  the  matrix 
will  be  sparse,  the  elements  cannot  be  arranged  in  such  a  way 
to  confine  the  elements  to  a  narrow  region  [10].  Although  the 
computation  time  will  depend  on  the  speed  of  the  computer,  in 
general  the  technique  may  be  too  time  consuming  when  dealing 
with  problems  that  require  high  resolution  and  a  large  number 
of  propagation  steps.  The  alternate  direction  method  was  devel¬ 
oped  to  address  this  difficulty  [9] . 


V.  Alternate  Direction  Implicit  Method 

The  alternate  direction  implicit  method  is  a  modification 
of  the  Crank-Nicolson  method  that  uses  smaller  tridiagonal 
matrices.  The  ADI  method  can  be  derived  directly  from  the 
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Fig.  2.  The  marching  planes  of  the  three-dimensional  Crank-Nicolson  method. 
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Fig.  3.  The  line  by  line  decomposition  of  the  planes  using  the  ADI  method. 
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where  0(*)  represents  the  order  of  error.  The  error  introduced 
by  the  new  error  term  is  0(Az 3)  and  can  be  ignored  because  it 
is  of  the  same  order  as  the  truncation  error.  Equation  (20)  above 
can  be  re-written  in  factor  form 


( 


i  - 


un+ 1 


(21) 


We  can  split  (21)  in  a  way  that  allows  us  to  solve  for  the  un¬ 
known  field  line  by  line  (direction  by  direction) 

(1-++’*I/2  =  (1+i+*)"'>  <22) 

H a7‘”Mi+4§7“"+1/2'  <23) 


Combined,  (22)  and  (23)  are  known  as  the  Peaceman-Rachford 
method  [8],  [9],  [17]  and  is  equivalent  to  (21)  (easily  verified 
by  multiplying  both  sides  by  (1  —  (rx/4jkG)8x)).  Fig.  3  shows 
schematics  of  the  line  by  line  decomposition  of  the  planes  within 
one  propagation  step.  For  this  example,  the  unknowns  of  the 
intermediate  plane  can  be  solved  using  one  Nx  matrix  and  the 
unknowns  of  the  n  + 1  plane  can  be  solved  using  one  Ny  matrix. 

One  difficulty  with  the  ADI  technique,  is  that  it  has  special 
conditions  for  boundary  values  in  the  n  +  1/2  plane.  The  n  and 
the  n+ 1  planes  are  physical  planes  and  the  intermediate  n+ 1/2 


plane  may  be  considered  a  virtual  plane  which  may  not  have  the 
same  boundary  conditions  as  the  physical  plane.  Adding  (22) 
and  (23),  we  arrive  at  an  expression  for  the  intermediate  plane 
in  terms  of  the  physical  planes 


2  V  4jka  y 


.71+1 
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u  =  f3(x,  y,  z)  =  boundary. 


»y)un  (24) 
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Equation  (24)  can  only  be  satisfied  if  the  boundary  values 
are  known  for  each  propagation  step.  If  we  impose  the  same 
boundary  conditions  on  the  intermediate  plane  as  the  physical 
plane,  then  u”+,1/2  =  fl£t1/2  =  =  ft*1,  and  the  overall 

accuracy  will  be  first  order  [9].  As  with  the  Crank-Nicolson 
method,  the  ADI  method  is  also  unconditionally  stable  [9], 
[10],  [17],  [29]. 

Tunnels  with  rectangular  cross  sections  are  the  simplest  for 
ADI  to  handle  because  one  matrix  can  be  used  to  solve  for  the 
field  along  the  x  axis  and  one  matrix  for  the  field  along  the  y 
axis.  However,  for  a  circular  case,  there  will  need  to  be  (Ny  —  2) 
matrices  to  solve  for  the  field  along  the  x  axis  and  (Nx  —  2) 
matrices  to  solve  for  the  field  along  the  y  axis.  The  Taylor  series 
approximation  is  used  to  approximate  the  curved  boundary  [10] 
and  different  ADI  matrices  must  be  used  for  each  line.  The  size 
of  the  matrices  will  depend  on  the  number  of  mesh  points  along 
a  line  of  the  circle. 

Another  complication  arises  when  enforcing  Neumann-type 
boundary  conditions.  Using  the  technique  outlined  in  Morton 
[10],  we  can  derive  a  first  order  normal  derivative  using  interpo¬ 
lation  in  order  to  express  boundary  mesh  points  in  terms  of  inte¬ 
rior  mesh  points.  As  a  result,  the  Neumann  boundary  condition 
may  couple  two  lines  in  one  dimension  and  increase  our  number 
of  unknowns.  This  makes  the  order  in  which  the  line  by  line  de¬ 
composition  occurs  important.  For  a  circle,  the  line  by  line  de¬ 
composition  must  begin  at  the  center  line  where  the  boundary 
values  are  reflected  inward  and  there  is  no  coupling.  This  makes 
the  complete  solution  of  the  center  line  possible  without  the  in¬ 
troduction  of  new  unknowns.  The  lines  directly  adjacent  to  the 
center  line  must  be  solved  next  and  so  on  until  the  field  along 
the  plane  is  solved  completely.  In  order  to  avoid  increasing  the 
number  of  unknowns  and  the  size  of  the  ADI  matrices,  the  lines 
must  be  solved  in  successive  order. 
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TABLE  I 

The  Normalized  RMS  Error  for  the  Rectangular,  Circular  and  Semi-Circular 
Cylindrical  Waveguides  With  Dirichlet  and  Neumann  Boundary  Conditions 


Erms  Error  (%) 

Rectangular  WG 

Circular  WG 

Semi-Circular  WG 

Ax,A  y 

Dirichlet 

Neumann 

Dirichlet 

Neumann 

Dirichlet 

Neumann 

0.8A 

5  A 

11.9 

10.3 

11.4 

25.2 

12.1 

56.4 

0.8A 

10  A 

13.7 

12.2 

13.2 

24.2 

14 

55.1 

0.8A 

20  A 

19.5 

18 

18.7 

21.9 

19.7 

51 

0.4A 

5  A 

4.9 

4.7 

5.0 

12.9 

5.2 

29.1 

0.4A 

10  A 

7.3 

7.1 

7.0 

11.9 

7.6 

27.5 

0.4A 

20  A 

14.4 

14.2 

13.8 

12 

14.6 

24.3 

VI.  PEC  Waveguide  Examples  Using  the  ADI  Method 

In  this  section,  the  analytical  solutions  of  simple  PEC  waveg¬ 
uides  are  used  to  measure  the  accuracy  of  the  ADI  approach. 
Our  analysis  will  consider  three  test  cases:  the  (a)  square,  (b) 
circular  and  (c)  semi-circular  cylindrical  PEC  waveguides. 

All  simulations  are  done  using  a  unit  strength  Gaussian  with 
3. 5 A  standard  deviation  placed  at  the  center  of  the  transverse 
plane  as  the  initial  field  and  operating  at  a  frequency  of  3  GHz. 
The  magnetic  and  electric  scalar  potentials  are  simulated  by  en¬ 
forcing  either  the  Dirichlet  or  Neumann  boundary  conditions. 
In  test  case  (a),  the  waveguide  cross-section  is  40 A  x  40A  and 
in  test  cases  (b)  and  (c)  the  radius  is  20 A.  All  test  cases  are  sim¬ 
ulated  up  to  a  distance  of  100  m. 

Under  these  parameters,  only  modes  propagating  within 
=bl5°  make  significant  contributions  and,  therefore,  the  PE  ap¬ 
proximation  is  valid.  Furthermore,  because  of  the  PE  limitation, 
Nyquist’s  theorem  restricts  the  descretizations,  Ax  and  Ay,  to 
be  less  than  1.9A.  This  is  shown  by  considering  the  transverse 
wavelength  at  15°  [7]. 

Table  I  displays  the  rms  error 
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Fig.  4.  The  (a)  analytical  solution  and  (b)  numerical  approximation  of  the  rect¬ 
angular  waveguide  with  Dirichlet  boundary  conditions  (RMS  error  =  7.3%) 
and  the  (c)  analytical  solution  and  (d)  numerical  approximation  of  the  rect¬ 
angular  waveguide  with  first  order  Neumann  boundary  conditions  (RMS  error 
=  7.1%)  . 
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where  is  the  approximated  discretized  field,  u^f  is  the 

known  analytical  field  and  N  is  the  total  number  of  unknowns. 
The  analytical  solutions  used  in  Sections  VI-A-C  are  obtained 
in  the  traditional  ways  outlined  by  Harringtion  [1]  and  other 
texts. 

A.  Test  Case  (a):  Rectangular  Waveguide  With  Dirichlet  and 
Neumann  Boundary  Conditions 

Fig.  4(a)-(b)  show  the  analytical  and  numerical  solutions  for 
the  rectangular  waveguide  with  Dirichlet  boundary  conditions. 
Fig.  4(c)-(d)  show  the  analytical  and  numerical  solutions  for  the 
rectangular  waveguide  with  Neumann  boundary  conditions.  It 
is  clear  from  these  figures  that  there  is  very  good  agreement  be¬ 
tween  the  analytical  and  simulated  field  patterns.  Table  I  shows 
that  for  transverse  descretizations,  Ax  =  Ay  —  0.4A,  and 
marching  step  size,  A z  =  5 A,  the  rms  error  is  within  5%. 


B.  Test  Case  (b):  Circular  Cylindrical  Waveguide  With 
Dirichlet  and  Neumann  Boundary  Conditions 

Fig.  5(a)-(b)  show  the  analytical  and  numerical  solutions 
for  the  circular  waveguide  with  Dirichlet  boundary  conditions. 
Again,  there  is  good  agreement  between  the  analytical  and 
numerical  field  patterns.  Table  I  shows  that  for  transverse 
descretizations,  Ax  —  Ay  =  0.4A,  and  marching  step  size, 
A z  =  5 A,  the  rms  error  is  within  5%.  Fig.  5(c)-(d)  show  the 
analytical  and  numerical  solutions  for  the  circular  waveguide 
with  Neumann  boundary  conditions,  using  first  order  interior 
interpolation.  There  is  good  agreement  between  the  analytical 
and  numerical  field  patterns  but  the  rms  error  increases  to 
12.9%.  The  rms  error  for  the  Neumann  case  is  greater  than  the 
Dirichlet  case  because  both  the  boundary  conditions  and  the 
interpolations  are  first  order. 

C.  Test  Case  (c):  Semi-Circular  Cylindrical  Waveguide  With 
Dirichlet  and  Neumann  Boundary  Conditions 

Fig.  6(a)-(b)  show  the  analytical  and  numerical  solutions 
for  the  semi-circular  waveguide  with  Dirichlet  boundary 
conditions.  As  with  the  previous  test  cases,  there  is  very 
good  agreement  between  the  analytical  and  simulated  field 
patterns.  Table  I  shows  that  for  transverse  descretizations, 
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(a)  (b) 


-2-1012  -2-1012 


(c)  (d) 


Fig.  5.  The  (a)  analytical  solution  and  (b)  numerical  approximation  of  the  cir¬ 
cular  waveguide  with  Dirichlet  boundary  conditions  (RMS  error  =  7.0%)  and 
the  (c)  analytical  solution  and  (d)  numerical  approximation  of  the  circular  wave¬ 
guide  with  Neumann  boundary  conditions  (RMS  error  =  11.90%). 


Fig.  6.  The  (a)  analytical  solution  and  (b)  numerical  approximation  of  the  semi¬ 
circular  waveguide  with  Dirichlet  boundary  conditions  (RMS  error  =  7.6%) 
and  the  (c)  analytical  solution  and  (d)  numerical  approximation  of  the  semi¬ 
circular  waveguide  with  Neumann  boundary  conditions  (RMS  error  =  27.5%). 


Ax  =  Ay  =  0.4A,  and  marching  step  size,  Az  =  5A,  the  rms 
error  is  about  5.2%. 

Fig.  6(c)-(d)  show  the  analytical  and  numerical  solutions  for 
the  semi-circular  waveguide  with  Neumann  boundary  condi¬ 
tions.  There  is  very  good  agreement  between  the  field  patterns 
in  both  cases.  However,  Table  I  shows  for  Ax  =  Ay  =  0.4 A 
and  Az  =  5A,  the  rms  error  is  very  large  (29%).  Fig.  7  shows 
the  analytical  (solid-line)  and  simulated  (dashed-line)  fields  at 
(jy  =  7t/2.  As  shown  in  the  figure,  the  error  is  concentrated  near 
y  —  0  and  the  field  pattern  is  preserved.  Like  in  the  circular 
Neumann  case,  the  rms  error  can  be  improved  by  decreasing 
the  mesh  descretizations.  Also,  we  can  see  from  Table  I  the  error 
grows  slightly  as  the  step  size  is  decreased  with  fixed  mesh  de¬ 
scretizations.  This  is  because  the  field  is  converging  towards  a 
field  not  exactly  matching  the  analytical  solution. 


Fig.  7.  The  magnitude  of  the  analytical  (blue,  solid)  and  simulated  (red, 
dashed)  field  at  x  —  0  for  PEC  tunnel  with  semi-circular  cross-section  and 
Neumann  boundary  conditions. 


Although  the  ADI  method  can  produce  larger  errors  than  the 
CN  method,  for  the  examples  considered,  the  fields  generated 
by  both  codes  are  nearly  identical.  Keeping  the  descretizations 
the  same,  the  error  between  the  two  methods  is  less  than  1%.  A 
discussion  on  the  accuracy  of  the  ADI  method  can  be  found  in 
[30]. 

VII.  Impedance  Waveguide  Examples 
Using  the  ADI  Method 

Tunnels  with  impedance  walls  can  be  very  difficult  to  solve 
and,  aside  from  some  special  cases,  closed  form  solutions  may 
not  exist  [6],  [14].  In  Sections  VII- A  and  B,  we  will  consider 
lossy  tunnels  with  known  analytical  solutions.  We  will  con¬ 
sider  two  such  cases,  the  rectangular  tunnel  with  small  surface 
impedance  and  the  circular  tunnel  with  a  linear  dipole  located 
at  its  center. 

A.  Rectangular  Tunnel 

Exact  analytical  solutions  for  the  fields  in  rectangular  tunnels 
with  lossy  walls  do  not  exist  [6],  [14].  In  [6],  a  lossy  tunnel  is 
treated  as  an  imperfect  metal  waveguide  with  the  relative  surface 
impedance  of  a  good  conductor 

z  =  zs  =  Rs  +  jXs  =  -A=.  (27) 

V  J&r 

If  the  magnitude  of  the  relative  surface  impedance  is  much  less 
than  unity,  then  perturbation  techniques  can  be  used  to  derive 
approximate  fields  [6].  Unlike  the  perfect  waveguide  case,  pure 
TE  and  TM  modes  will  not  always  propagate  within  the  im¬ 
perfect  waveguide.  Instead,  there  will  exist  hybrid  (ExH)m  n 
and  ( EHX )m  n  modes  that  are  coupled  TE  and  TM  modes.  The 
(. ExH)mn  mode  is  a  TE  mode  with  a  small  TM  component, 
and  the  ( EHX  )m  n  mode  is  a  TM  mode  with  a  small  TE  com¬ 
ponent. 

_  In  the  case  of  the  rectangular  tunnel,  the  impedance  matrix 
T0GT0 ,  shown  in  (15),  becomes  diagonal  because  the  normal 
vectors,  nx  and  ny,  disappear  alternatively  on  the  tunnel  wall. 
Therefore,  the  transverse  reduced  field  components  U (x,  y ,  z) 
and  V(x,y,z)  in  the  orthogonal  (x,y)  plane  are  decoupled 
and  can  be  solved  independently  [4].  The  hybrid  fields  that  are 
solved  numerically  will  depend  on  the  boundary  conditions 
enforced  on  the  waveguide  walls. 
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(a)  EHx  (b)  ExH 

dE 


Fig.  8.  Rectangular  waveguide  profiles  displaying  the  boundary  conditions  for 
the  Ex  field  of  the  EHX  mode  and  the  Ey  field  of  the  EXH  mode  in  the  trans¬ 
verse  plane. 


Consider  the  simple  example  of  a  rectangular  waveguide  0  < 
x  <  H,  0  <  y  <  W,  with  three  PEC  walls  and  one  impedance 
wall  located  at  x  =  H ,  where  the  normal  vectors  are  nx  —  1  and 
ny  —  0.  Fig.  8  shows  the  boundary  conditions  at  the  three  PEC 
walls  for  the  EHX  [Fig.  8(a)]  and  the  EXH  [Fig.  8(b)]  hybrid 
modes  with  the  Ex  and  Ey  fields  given  by  [6,  Eq.  (36)]. 

Eet  us  consider  the  case  of  the  EHxn  mode  propagating 
in  the  waveguide  at  a  frequency  of  1  GHz  with  a  wall  of 
impedance,  =  0.01  +  j 0.01.  Good  conductors  have  large 
conductivity,  and  for  this  case  the  equivalent  conductivity 
is  556  S/m.  Fig.  9(a)  shows  the  magnitude  of  the  analytical 
solution  of  Ex  in  the  origin  plane  of  z  =  0  and  Fig.  9(b) 
shows  the  magnitude  of  the  analytical  solution  of  Ex  at  the 
distance  of  z  =  1  km.  Fig.  9(c)  shows  the  magnitude  of  the 
simulated  Ex,  using  the  analytical  field  in  the  origin  plane 
as  the  initial  field.  The  discretizations  along  the  transverse 
plane  are  Ax  =  Ay  =  0.27A,  and  along  the  propagation 
axis,  A z  =  11  A.  We  can  see  that  there  is  strong  agreement 
between  the  analytical  field  and  the  simulated  field  in  the 
transverse  plane.  The  field  pattern  doesn’t  depend  on  distance 
because  a  single  mode  was  used  as  the  initial  field.  The  ADI 
simulation  produced  the  same  pattern  throughout  the  step  by 
step  simulation.  For  this  case,  this  difference  between  the  mode 
attenuation  factors  (MAFs)  of  the  analytical  and  numerical 
solution  is  2.36  dB/km. 

Now,  let  us  consider  the  case  of  the  ExHn  mode  propa¬ 
gating  in  the  waveguide  at  a  frequency  of  1  GHz  with  a  wall 
of  impedance,  =  0.3  +  j 0.3.  For  this  case  the  equivalent 
conductivity  is  0.62  S/m.  Fig.  9(d)  shows  the  magnitude  of  the 
analytical  solution  of  Ey  in  the  origin  plane  of  z  =  0.  Fig.  9(e) 
shows  the  magnitude  of  the  analytical  solution  of  Ey  at  the  axial 
distance  of  z  =  5  km.  The  distance  of  5  km  was  chosen  so  that 
there  would  be  noticeable  attenuation  along  the  axial  direction. 
Fig.  9(f)  shows  the  magnitude  of  the  simulated  Ey  field,  using 
the  analytical  field  in  the  origin  plane  as  the  initial  field.  The 
discretizations  along  the  transverse  plane  are  the  same  as  in  the 
EHX n  example,  Ax  =  Ay  =  0.27A  and  Az  =  11A.  Again, 
as  in  the  EHxn  example,  we  can  see  that  there  is  strong  agree¬ 
ment  between  the  analytical  and  the  simulated  field  patterns  in 
the  transverse  plane.  Also,  the  magnitudes  of  the  analytical  and 
simulated  fields  show  that  the  simulated  field  captures  the  atten¬ 
uation  of  the  field  along  the  jz  direction. 


(a):  EHX  Analytical  field  at  z=0m 


(b):  Analytical  field  at  z=1km 


0 
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(c):  Numerical  field  at  z=1km 


(d):  EXH  Analytical  field  at  z=0m 
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(e):  Analytical  field  at  z=5km 
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(f):  Numerical  field  at  z=5km 
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Fig.  9.  The  magnitudes  of  the  analytical  and  numerical  fields  of  the 
(a-c)  EHxll  and  (d-f)  EXHX1  mode. 


Fig.  10.  Cross-section  of  a  circular  tunnel  with  lossy  walls  and  a  linear  source 
located  at  the  center. 


The  main  advantage  of  the  ADI  method  over  the 
Crank-Nicolson  method  is  the  use  of  smaller  ADI  matrices.  In 
the  rectangular  case,  Ex  and  Ey  fields  are  decoupled  and  each 
field  can  be  treated  as  a  scalar  PE  problem.  As  with  the  previous 
scalar  PE  problems,  if  there  are  Nx  and  Ny  mesh  points  along 
the  x  and  y  dimensions,  respectively,  the  Crank-Nicolson 
matrix  will  be  a  sparse  matrix  of  rank  (NxNy  —  4).  However, 
the  ADI  method  will  only  require  2  matrices  of  rank  Nx  and 
Ny,  respectively. 

B.  Circular  Tunnel 

Next,  we  use  the  ADI-PE  to  solve  for  the  field  in  a  lossy  tunnel 
with  a  circular  cross-section  of  radius  a  (shown  in  Fig.  10).  The 
tunnel  is  treated  as  a  hole  in  a  lossy  dielectric  medium  where 
the  interior  is  free  space  and  the  exterior  is  characterized  by 
(jj,  o,  62),  with  62  =  60  (er  —  j  a  r).  Fig.  10  also  shows  the  linear 
dipole  current  source  with  dipole  moment  P  (in  amp-meters), 
Jx  =  P8{(j))8{z)8{p) / p,  located  at  the  center  of  the  tunnel 
(p  =  0, 0  =  0)  in  the  2  =  0  plane.  The  derivation  of  the  electric 
fields  produced  by  this  source  is  well  known  and  the  details  are 
given  in  [14]  and  [15]. 
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(a)  p=  0.26667a  (b)  p=  0.93333a 


Fig.  11.  (a)-(b)  Axial  distribution  of  the  analytical  (blue,  solid)  and  simulated 
(red,  dashed)  Ex  field  intensity  for  the  H E1^1_2  modes;  a  =  2  m,/  =  1  GHz, 
P/a  —  1  amp,  er  =  12,  o 0  —  0.02  S/m. 


As  with  the  imperfect  rectangular  waveguide,  the  lossy  cir¬ 
cular  tunnel  will  produce  hybrid  modes.  For  the  special  case  of 
the  centered  linear  source,  the  only  possible  hybrid  modes  are 
EHim  and  HEim  [14],  [15].  The  transverse  electric  field  com¬ 
ponents  are  shown  in  [15,  Eq.  (95)  and  (96)]. 

For  our  ADI  simulations,  we  consider  a  smooth,  straight 
tunnel  of  radius  2  m  and  1  km  in  length.  We  choose  our 
electrical  parameters  to  match  the  parameters  studied  in  the 
published  work  of  Dudley  and  Mahmoud  [15].  We  choose 
er  =  12  and  aQ  =  0.02  S/m  and  P/a  =  1  amp.  In  this  case, 
HEn  is  the  dominant  mode  and  attenuates  least  over  the 
axial  distance  of  the  guide.  The  MAF  of  the  HEu  mode  is 
2.77  dB/100  m  [14].  In  a  lossy  tunnel,  higher  order  modes  will 
attenuate  more  than  the  dominant  mode  and,  for  large  axial 
distances,  the  field  intensity  will  correspond  to  the  intensity  of 
the  HEn  mode. 

Fig.  ll(a)-(b)  shows  the  analytical  field  intensity  (in  dB)  as 
a  function  of  axial  distance  at  </>  =  7t/2  for  about  (a)  p  =  0.25a 
and  (b)  p  =  0.9a  (in  blue).  The  first  two  modes  (HEi^-2)  are 
plotted  because  the  angles  of  propagation,  with  respect  to  the  z 
axis,  are  below  the  parabolic  approximation  limit  of  ±15°.  The 
angles  of  the  modes  are  calculated  using 


cos  0  = 


Rejkz) 


(28) 


The  complex  kz  terms  are  calculated  using  the  procedure 
outlined  in  [15,  (Eq.  76)].  The  ordering  of  the  first  three 
modes,  according  to  angle,  is  HEn(3.3°),  HEi2(7.3°)  and 
HEis(11.8°).  We  can  see  from  Fig.  11  that  the  field  inten¬ 
sity  can  be  divided  into  a  near-zone  and  a  far-zone  [14].  The 
near-zone  is  characterized  by  rapid  variations  due  to  higher 
mode  interactions.  The  far-zone  is  smoother  and  more  linear 
because  the  higher  order  modes  disappear  due  to  attenuation. 
In  the  far-zone  the  slope  corresponds  to  the  MAF  of  the  HEu 
mode.  For  our  analysis  we  define  the  Erms  as  a  function  of 
axial  distance 


Erms  =  20  log 
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The  analytical  solution  for  the  HE11  and  HE\^  modes  at 
the  2  =  0  plane  is  used  as  the  initial  field  for  the  ADI  simula¬ 
tion  and  the  discretizations  along  the  transverse  plane  and  the 
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Fig.  12.  (a)-(b)  Analytical  transverse  Ex  and  Ey  field  distribution  at  a  distance 
of  1  km  for  the  HE1^1_2  modes,  respectively,  (c)-(d)  Numerical  transverse  Ex 
and  Ey  field  distribution  with  A x  —  0.444A,  Ay  —  0.444A,Az  =  1.6667A. 


axial  direction  are  Ax  =  Ay  =  0.444A  and  Az  =  1.667A. 
As  shown  in  Fig.  11,  the  analytical  and  simulated  results  are  in 
very  good  agreement.  In  the  far-zone  (z  >  600  m),  where  the 
field  becomes  more  linear  the  Erms  doesn’t  exceed  —10  dB. 
As  we  can  see  from  Fig.  11,  the  trend  and  shape  of  the  analyt¬ 
ical  result  is  captured  in  the  simulation.  Fig.  12  compares  the 
magnitudes  of  the  Ex  and  Ey  field  distribution  in  the  transverse 
plane.  We  can  see  that  the  magnitudes  and  field  patterns  are  in 
close  agreement  between  the  analytical  [Fig.  12(a)-(b)]  and  nu¬ 
merical  [Fig.  12(c)-(d)]  results. 

Unlike  the  rectangular  case,  the  impedance  matrix  T0G0T0 
(15),  is  not  diagonal  and  the  Ex  and  Ey  fields  have  to  be  solved 
simultaneously.  As  a  result,  we  have  to  solve  a  true  vector  PE 
and  the  size  of  both  the  Crank-Nicolson  and  ADI  matrices  in¬ 
crease.  If  there  are  Nx  and  Ny  mesh  points  along  the  x  and  y 
dimensions,  respectively,  the  Crank-Nicolson  matrix  will  be  a 
sparse  matrix  of  rank  2  x  (NxNy  —  4),  double  the  size  of  the 
scalar  PE  case.  The  ADI  method  will  require  Ny  —  2  matrices 
that  are  at  most  size  2  xNx  and  Nx  —  2  matrices  that  are  at  most 
size  2  xNy.  When  dealing  with  large  problems,  the  memory  re¬ 
quired  to  perform  calculations  with  the  Crank-Nicolson  matrix, 
even  with  use  of  the  sparse  command  in  MATFAB,  might  be¬ 
come  too  computationally  intensive  for  an  average  PC. 

VIII.  Comparing  ADI  Simulations  With 
Experimental  Data 

Finally,  we  consider  the  case  of  real  tunnels  and  compare 
ADI  simulations  to  measurements  shown  in  published  results. 
Our  first  example  is  the  Massif  Central  tunnel  in  south-central 
France  studied  in  the  published  work  of  Dudley,  Fienard,  Mah¬ 
moud,  and  Degauque  [14].  The  tunnel  is  straight  and  3.5  km 
in  length  with  smooth  walls  composed  of  large  blocks  of  stone 
or  concrete.  The  roughness  is  estimated  to  be  in  the  order  of  a 
centimeter.  The  transverse  dimensions  of  the  tunnel  are  given  in 
Fig.  13.  The  experiment  performed  by  Dudley  et  al.  consisted 
of  a  transmitting  and  receiving  antenna  placed  vertically  at  a 
height  of  2  m  and  horizontally  at  one-quarter  the  width  of  the 
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Fig.  13.  The  profile  of  the  Massif  Central  tunnel. 


tunnel  (z  =  2.15  m),  shown  as  an  “X”  in  Fig.  13.  The  trans¬ 
mitting  power  was  34  dBm  and  the  frequencies  of  450  MHz 
and  900  MHz  are  studied.  Half-wave  vertical  dipoles  are  used 
for  450  MHz, and  vertically  polarized  wide-band  horn  antennas 
(with  a  gain  of  7  dBi)  are  used  for  900  MHz. 

The  research  team  of  Dudley  et  al.  [14]  uses  an  equivalent 
rectangular  tunnel  (see  Fig.  13)  in  order  to  characterize  the 
Massif  Central  tunnel.  The  dimensions  of  the  rectangular  tunnel 
are  shown  in  Fig.  13  and  the  electrical  parameters  are  chosen 
to  be  er  =  5  and  aQ  =  0.01  S/m.  The  theoretical  attenuation 
constant  of  the  dominant  mode  of  the  equivalent  rectangular 
tunnel  is  given  by  [13]  and  [14] 

M2 


OL 11  = 


H  \2H)  Re  Lv^ 


—  ( -  V 

+  W  \2WJ 


Re 


y/Cpc  1 
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where  H  is  the  horizontal  length  and  W  is  the  vertical  length, 
given  by  H  =  7.8  m  and  W  =  5.3  m.  The  attenuation  of 
the  dominate  mode  in  the  equivalent  rectangular  tunnel  corre¬ 
sponds  to  the  attenuation  of  the  measured  field  in  the  Massif 
Central  tunnel.  Figs.  14  and  15  show  the  field  intensity  in  dBm 
as  a  function  of  axial  distance  (in  blue)  as  measured  by  the  re¬ 
search  team  of  Dudley  et  al.  [14].  The  MAF  of  the  equivalent 
rectangular  tunnel  and  the  experimental  data  are  shown  in  the 
first  two  rows  of  Table  II. 

The  ADI  simulations  are  done  using  the  actual  arched  cross- 
section  of  the  Massif  Central  tunnel  and  the  equivalent  rectan¬ 
gular  tunnel  (Fig.  13).  The  electrical  parameters  are  chosen  to 
match  the  parameters  of  the  research  team.  At  this  point,  it  is 
important  to  discuss  possible  sources  of  error  in  the  simulation 
of  real  tunnels.  One  possible  error  may  be  the  incorrect  predic¬ 
tion  of  electrical  parameters  of  the  tunnel  walls.  If  a  real  tunnel 
is  made  up  of  dry  cement  and  stone,  the  parameters  will  depend 
on  the  materials  used  to  mix  the  cement  and  its  percentage  of 
water  content.  One  source  estimates  the  parameters  of  dry  con¬ 
crete  to  vary  from  er  =  4  —  10  and  aQ  =  0.001  —  0.01  S/m 
and  for  wet  concrete,  er  =  10  —  20  and  a0  =  0.01  —  0.1  S/m 
[16].  In  addition,  the  electrical  parameters  will  also  change  with 
frequency.  For  the  equivalent  rectangle,  this  isn’t  a  problem  be¬ 
cause  the  parameters  have  been  fixed  in  [14]. 

The  initial  field  used  in  the  simulations  may  not  represent  the 
actual  antenna  field  used  in  measurements.  However,  even  with 
an  incorrect  initial  field,  the  field  in  the  far  zone,  where  lower 
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Fig.  14.  The  electric-field  intensity  in  dBm  and  the  corresponding  least  square 
fit  line  of  the  (a)  Massif  Central  tunnel  (blue,  solid),  (b)  the  ADI  solution  using 
the  equivalent  rectangular  profile  with  er  —  5  and  o0  —  0.01  S/m  (red, 
dashed),  (c)  the  ADI  solution  using  the  arched  profile  with  er  =  5  and  cr0  — 
0.01  S/m  (green,  circle)  and  (d)  the  ADI  solution  using  the  arched  profile  with 
er  —  12  and  o0  —  0.01  S/m  (black,  asterisk)  as  a  function  of  axial  distance 
and  operating  frequency  of  450  MHz. 


Fig.  15.  The  electric-field  intensity  in  dBm  and  the  corresponding  least  square 
fit  line  of  the  (a)  Massif  Central  tunnel  (blue,  solid),  (b)  the  ADI  solution  using 
the  equivalent  rectangular  profile  with  er  =  5  and  aQ  —  0.01  S/m  (red, 
dashed),  and  (c)  the  ADI  solution  using  the  arched  profile  with  er  —  5  and 
<7 0  —  0.01  S/m  (green,  circle)  as  a  function  of  axial  distance  and  operating 
frequency  of  900  MHz. 


order  modes  dominate,  the  correct  attenuation  will  be  obtained. 
This  allows  us  some  flexibility  when  setting  up  our  simulations. 
So  long  as  the  lower  order  modes  are  illuminated,  different  ini¬ 
tial  fields  will  yield  the  same  MAFs  in  the  far  zone.  This  dis¬ 
crepancy  is  accounted  for  by  offsetting  the  vertical  axis  of  the 
simulation  data  so  the  least  square  fitted  lines  intersect  at  the 
far  zone  point.  Instead  of  having  the  dipole  inside  the  tunnel  en¬ 
trance,  as  in  the  Massif  Central  experiment,  we  set  up  the  dipole 
30  m  away  from  the  tunnel  opening.  This  allows  us  to  use  the 
far  field  expressions  for  the  dipole  and  avoid  solving  more  com¬ 
plicated  field  patterns.  We  do  not  know  the  value  of  the  field 
along  the  boundary  wall  in  the  initial  plane.  If  we  use  the  free 
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TABLE  II 

The  Experimental,  Theoretical  and  Simulated  Mode  Attenuation 
Factors  for  the  Massif  Central  Tunnel 


Frequency  (MHz)  /  Far  zone  range 

450  (.25-2. 5km) 

900(0.5-2. 5km) 

MAF  data(dB/km) 

33.0 

8.5 

■MAF  equiv  .rectangle  (dB  /  klTlj 

35.1 

8.6 

MAF AD j Arch  (dB/km) 

24 

8.45 

MAF AD  I  rectangle  (dB/km) 

35.3 

9.1 

space  values  along  the  boundary,  the  error  will  be  in  the  form  of 
high  oscillations  in  the  MAFs  in  the  far  zone.  We  can  eliminate 
these  oscillations  by  filtering  off  the  field  at  the  boundary  with 
a  unit  gaussian  with  a  =  0.75  m  (450  MHz)  and  a  —  1.5  m 
(900  MHz).  We  don’t  expect  an  exact  match  between  the  simu¬ 
lated  results  and  the  measured  results  in  the  near  zone  where  the 
parabolic  approximation  is  not  satisfied.  However,  we  expect,  in 
the  far-zone,  the  E-field  and  the  field  attenuation  to  be  correctly 
modeled  by  the  simulation.  Figs.  14-15  show  the  simulated  field 
intensity  in  dBm  of  the  tunnel  with  an  arched  cross  section  (in 
green,  circle)  and  the  equivalent  rectangle  (in  red,  dashed).  The 
MAFs  of  the  simulations  are  computed  from  the  slope  of  the 
least  square  fit  line  using  data  in  the  far  zone.  The  MAFs  of  the 
simulations  and  the  far  zone  ranges  are  shown  in  Table  II. 

In  Fig.  14,  for  the  450  MHz  case,  we  can  see  that  there  is 
good  agreement  between  the  simulated  rectangular  tunnel  and 
the  measured  data.  Table  II  shows  the  discrepancies  of  the  MAFs 
of  the  measured  data,  the  equivalent  rectangle  and  the  simu¬ 
lated  rectangle  are  about  2  dB/km.  The  discrepancy  of  the  sim¬ 
ulated  results  for  the  arched  shape  tunnel  and  the  measured  data 
is  about  11  dB/km.  This  is  believed  to  be  the  result  of  using 
an  incorrect  electrical  parameter.  Increasing  the  magnitude  of 
the  electrical  parameters  for  the  arched  shape  will  increase  the 
attenuation  of  the  tunnel.  This  is  shown  in  Fig.  14,  where  the 
black  curve  is  the  simulation  for  er  =  12  and  cr0  =  0.01  S/m. 
In  Fig.  15,  for  the  900  MHz  case,  we  see  that  there  is  more 
agreement  among  all  the  simulations.  As  Table  II  shows,  the  dis¬ 
crepancies  among  all  the  simulations  are  within  1  dB/km.  The 
good  agreement  may  be  due  to  the  frequency  dependant  nature 
of  the  electrical  parameters.  At  900  MHz,  the  electrical  parame¬ 
ters  used  for  the  equivalent  rectangle  may  be  in  better  agreement 
with  the  actual  electrical  parameters  of  the  tunnel  wall. 

Our  second  example  is  a  railway  tunnel  in  Japan  studied  in  the 
published  work  of  Chiba,  Inaba,  Kuwamoto,  Banno,  and  Sato 
[12].  The  tunnel  is  straight  and  1.47  km  in  length  and  also  made 
up  of  concrete.  The  cross-section  of  the  tunnel  is  arched  shaped 
and  has  two  3.9  m  lanes  separated  by  a  notch.  Fig.  16  shows 
the  dimensions  of  the  tunnel  cross-section.  The  measurement 
team  used  half-wave  dipole  antennas  for  the  transmitters  and  re¬ 
ceivers  and  the  antennas  were  placed  at  a  height  of  3.6  m  above 
the  ground  and  30  m  from  the  tunnel  opening  with  a  transmitting 


(a)  (b) 

Fig.  16.  The  profiles  of  the  (a)  Japanese  National  Railway  tunnel  and  (b)  the 
equivalent  circular  tunnel. 


Fig.  17.  The  (a)  MAF  of  the  Japanese  National  Railway  tunnel  (blue,  solid),  (b) 
the  electric- field  intensity  and  least  square  fit  line  of  the  ADI  solution  using  the 
equivalent  circular  profile  (red,  asterisk)  and  (c)  the  electric-field  intensity  and 
least  square  fit  line  of  the  ADI  solution  using  the  arched  profile  (green,  circle) 
as  a  function  of  axial  distance  with  an  operating  frequency  of  470  MHz. 

power  of  1  W.  The  measurement  team  characterized  the  tunnel 
with  an  equivalent  area  circular  cylindrical  waveguide  of  radius 
4.2  m.  The  electrical  parameters  of  the  equivalent  circular  cylin¬ 
drical  tunnel  were  chosen  so  that  the  attenuation  of  the  dominant 
mode  will  correspond  to  the  attenuation  of  the  measured  field. 
The  theoretical  attenuation  constant  is  given  by  (31),  shown  at 
the  bottom  of  the  page,  where  unrn  is  the  root  of  the  (n  —  l)th 
Bessel  function.  The  MAF  of  the  equivalent  circular  cylindrical 
tunnel  is  defined,  er  =  5.5  and  a0  =  0.03  S/m.  Figs.  17  and  18 
show  the  theoretical  MAF  of  the  dominant  mode  of  the  equiva¬ 
lent  circular  cylindrical  tunnel  (blue). 

The  ADI  simulations  are  done  using  the  arched  shape  cross 
section  (Fig.  16)  and  the  equivalent  circular  cross  section.  As 
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Fig.  18.  The  (a)  MAF  of  the  Japanese  National  Railway  tunnel  (blue,  solid), 
(b)  the  electric-field  intensity  and  least  square  fit  line  of  the  ADI  solution  using 
the  equivalent  circular  profile  (red,  asterisk)  and  (c)  the  electric-field  intensity 
and  least  square  fit  line  of  the  ADI  solution  using  the  arched  profile  (green, 
circle)  as  a  function  of  axial  distance  with  an  operating  frequency  of  900  MHz. 


TABLE  III 

The  Experimental,  Theoretical  and  Simulated  Mode 
Attenuation  Factors  for  the  Railway  Tunnel 


Frequency  (MHz) 

470  (.25-2. 5km) 

900(0. 5-5km) 

MAF  data(dB/km) 

9.8 

2.0 

MAF ADIequivalentcircle  {dB  j km) 

10.8 

2.95 

MA¥ADIArch  {dB/ km) 

16.32 

5.3 

before,  the  electrical  parameters  are  chosen  to  match  the  pa¬ 
rameters  of  the  research  team.  Figs.  17  and  18  show  the  sim¬ 
ulated  field  intensity  of  the  tunnel  with  an  arched  cross-section 
(in  green,  circle)  and  the  equivalent  circle  (in  red,  asterisk).  The 
field  intensity  is  shown  in  dB  such  that  1  /TV/ m  =  0  dB. 

In  Fig.  17,  for  the  470  MHz  case,  we  can  see  that  there  is 
reasonable  agreement  between  the  attenuation  of  the  simulated 
equivalent  circular  cylindrical  tunnel  and  the  theoretical  atten¬ 
uation.  As  with  the  previous  tunnel,  the  simulation  is  least  ac¬ 
curate  in  this  frequency  range.  As  we  increase  the  frequency 
there  is  greater  agreement  between  the  theoretical  MAF  and  the 
simulated  fields.  As  before,  we  can  see  from  the  figures  that  in 
the  470  MHz  and  900  MHz  cases,  the  two  mode  interactions  are 
present  in  the  far  zone  [14].  Table  III  summarizes  the  theoretical 
MAF,  the  simulated  equivalent  circular  cylindrical  tunnel  and 
the  tunnel  with  an  arched  cross  section  MAFs.  The  discrepancy 
between  the  theoretical  MAF  and  the  equivalent  circular  cylin¬ 
drical  tunnel  is  about  1  dB/km  and  lower  for  both  frequencies. 
The  discrepancy  between  the  theoretical  MAF  and  the  tunnel 
with  arched  cross  section  about  6.5  dB/km  for  470  MHz  and 
3.3  dB/km  for  900  MHz. 

A.  Usage  of  CPU  Time  and  Memory 

In  this  section,  we  briefly  examine  the  computation  time 
and  memory  usage  of  the  Crank-Nicolson  and  ADI  method. 
Table  IV  shows  a  side  by  side  comparison  of  the  execution  time 
and  CPU  memory  of  the  Crank-Nicolson  and  ADI  method 


TABLE  IV 

The  CPU  Memory  and  Computational  Time  for  Nx  —  70  and 
2000  Marching  Steps  for  the  Tunnel  With  Rectangular 
Profile  and  Dirichlet  Boundary  Conditions 


Rectangle  Dirichlet 

ADI 

CN 

Execution  Time  (sec) 

124  (2min) 

3532  (lhr) 

Memory 

3.3KB  (208  elements  ) 

378.6KB  (23662  elements) 

for  the  case  of  a  PEC  tunnel  with  a  rectangular  profile  with 
resolution  Nx  =  Ny  =  70  and  using  2000  marching  steps. 
The  execution  time  refers  to  the  time  it  take  MATLAB  to 
perform  the  one  Backslash  command  and  multiplied  by  the 
number  of  loops.  This  accounts  for  any  variation  in  codes  by 
only  recording  the  time  it  takes  to  invert  the  matrices.  For 
the  Crank-Nicolson  method,  the  Backslash  command  was 
performed  after  the  memory  saving  Sparse  command  was  used. 
As  we  can  see  from  Table  IV,  ADI  runs  mush  faster  and  takes 
significantly  less  memory  than  the  Crank-Nicolson  method. 

To  continue  the  comparison  between  the  Crank-Nicolson 
method  and  the  ADI,  we  will  breifly  address  the  issue  of 
stability.  Following  the  discussion  in  [29],  on  the  stability  of 
the  ADI,  solving  for  a  step/wavenumber  independent  bound  on 
the  amplification  matrix  becomes  difficult  because  the  ADI-PE 
matrices  contain  off-diagonal  terms  that  make  finding  eigen¬ 
values  difficult.  For  the  PEC  waveguide,  numerical  tests  have 
shown  that  the  ADI  amplification  matrix  is  normal  and  bounded 
by  unity  and  the  von  Neumann  condition  is  sufficient.  When 
the  impedance  boundary  condition  is  applied,  the  matrices  are 
no  longer  normal,  but  the  amplification  matrix  is  still  bounded 
by  unity  for  various  discretizations  and  steps.  For  examples 
considered  in  the  paper,  we  have  not  encountered  instability 
when  applying  the  IBC  to  the  ADI-PE.  The  CN  is  uncondition¬ 
ally  stable  and  less  sensitive  to  boundary  conditions.  However, 
for  modest  discretizations,  we  were  able  to  get  very  good 
agreement  between  the  CN  and  ADI  solutions,  which  would 
imply  that  the  ADI  is  also  stable  with  the  IBC. 

IX.  Conclusion 

We  have  presented  the  ADI  method  for  use  in  solving  the 
parabolic  equation  for  radiowave  propagtion  in  tunnels.  We  have 
shown  that  because  the  ADI  technique  uses  smaller  matrices 
than  the  widely  used  Crank-Nicolson  method,  it  reduces  the 
computational  labor  of  an  average  PC. 

Using  analytical  models,  we  have  shown  that  the  ADI 
method,  in  combination  with  the  parabolic  approximation,  can 
simulate  the  attenuation  and  field  patterns  in  tunnels  whose 
electrical  parameters  are  known.  Our  simulation  results  also 
verify  the  use  of  the  grazing  angle  impedance  approximation 
on  the  tunnel  walls. 

In  the  case  of  actual  tunnels,  our  ADI-PE  approach  com¬ 
pares  well  with  published  experimental  data  for  the  900  MHz 
frequency  in  the  far  zone  for  the  Massif  Central  tunnel.  The 
tunnel’ s  far  zone  is  the  region  where  rapid  oscillations  cease  and 
the  lowest  order  modes  dominate.  The  descrepancies  in  mea¬ 
sured  and  simulated  results  may  be  attributed  to  lack  of  knowl¬ 
edge  of  the  tunnel’s  electrical  parameters.  Simulated  ADI-PE 
results  for  equivalent  tunnel  profiles,  where  electrical  parame¬ 
ters  are  defined,  show  high  accuracy. 
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Future  work  will  use  the  ADI-PE  to  study  the  cases  of  curved 
and  branching  tunnels  with  rough  surfaces  as  well  as  making 
comparisons  to  experimental  data. 
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The  four-state  random  walk  (4RW)  model,  wherein  the  particle  is  endowed  with 
two  states  of  spin  and  two  states  of  directional  motion  in  each  space  coordinate, 
permits  a  stochastic  solution  of  the  Schrodinger  equation  (or  the  equivalent  para¬ 
bolic  equation)  without  resorting  to  the  usual  analytical  continuation  in  complex 
space  of  the  particle  trajectories.  Analytical  expressions  are  derived  here  for  the 
various  transitional  probabilities  in  a  4RW  by  employing  generating  functions  and 
eigenfunction  expansions  when  the  particle  moves  on  a  1  + 1  space-time  lattice  with 
two-point  boundary  conditions.  The  most  general  case  of  dissimilar  boundaries 
with  partially  reflecting  boundary  conditions  is  treated  in  this  paper.  The  transi¬ 
tional  probabilities  are  all  expressed  in  terms  of  a  finite  summation  involving  trigo¬ 
nometric  functions  and/or  Chebyshev  polynomials  of  the  second  kind  that  are  char¬ 
acteristics  of  diffusion  and  Schrodinger  equations,  respectively,  in  the  4RW  model. 
Results  for  the  special  case  of  perfectly  absorbing  boundaries  are  compared  to 
numerical  values  obtained  by  directly  counting  paths  in  the  random  walk 
simulations.  ©  2009  American  Institute  of  Physics.  [DOI:  10.1063/1.3122768] 


I.  INTRODUCTION 

The  four-state  random  walk  (4RW)  model,  wherein  a  particle  undergoing  random  walk  is 
endowed  with  two  states  of  direction  (in  a  one-dimensional  case)  and  two  states  of  spin  or  parity, 
was  considered  by  Ord  and  Deakin1  to  arrive  at  a  macroscopic  model  for  the  Schrodinger  equation 
and  physically  interpret  its  wavelike  solutions.  It  was  shown  in  that  paper  that  both  the  traditional 
diffusion  equation  as  well  as  the  Schrodinger  equation  were  embedded  in  the  same  physical 
model.  The  usual  diffusion  process  is  contained  in  the  overall  sum  of  all  particle  paths  irrespective 
of  their  direction  and  parity,  while  the  wavelike  behavior  of  the  Schrodinger  equation  is  contained 
in  the  differences  in  densities  of  the  right-going  particles  or  left-going  particles  with  opposite 
parity.  The  4RW  model  owes  its  existence  to  the  Feynman  chessboard  model  as  elaborated  in  Ref. 
2  and  is  also  useful  in  solving  practical  electromagnetic,  acoustic,  and  optical  boundary  value 
problems  for  the  complex  field  amplitude  when  a  stochastic  approach  is  used  to  treat  the  govern¬ 
ing  parabolic  equation.  (In  the  applied  sciences  area,  the  parabolic  equation  is  sometimes  referred 
to  as  the  parabolic  wave  equation.)  The  standard  parabolic  equation  used  in  such  time-harmonic 
problems  contains  partial  derivatives  with  respect  to  the  spatial  coordinates  only,  where  the  spatial 
coordinate  along  the  axial  direction  takes  the  place  of  the  time  variable  in  the  time-dependent 
Schrodinger  equation.  Furthermore,  the  potential  function  present  in  the  Schrodinger  equation  is 
replaced  by  the  medium  refractive-index  term  in  the  parabolic  equation.  The  parabolic  equation  is 
obtained  when  the  Helmholtz  equation  describing  the  true  field  is  subject  to  a  one-way  propaga¬ 
tion  with  a  subsequent  application  of  the  paraxial  approximation.3,4  Normally  one  needs  to  resort 
to  analytical  continuation  of  boundary  data,  as  is  done  in  Ref.  5,  when  the  parabolic  equation  is 
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solved  numerically  using  a  stochastic  approach.  The  resulting  random  walks  will  then  traverse  a 
complex-valued  space,  which,  in  turn,  calls  for  analytical  continuation  of  boundary  data  and  the 
spatial  geometry.6  However,  the  4RW  model  permits  a  solution  to  these  problems  without  such 
analytical  continuations. 

In  a  previous  paper,7  the  author  developed  expressions  for  the  various  transitional  probabilities 
for  the  4RW  on  a  discrete  lattice  for  spatially  unbounded  case.  In  addition  to  detailing  a  method 
for  arriving  at  the  transitional  probabilities,  the  results  presented  therein  could  also  be  useful  in 
determining  other  stochastic  quantities  of  interest  such  as  the  first  passage  time  probabilities, 
expected  number  of  visits  to  a  given  site,  and  maximum  excursions  of  a  random  walk  on  a  line  for 
various  linear  combinations  of  these  probabilities.  In  developing  numerical  schemes  for  solving 
boundary  value  problems  with  complex  geometries,  it  is  desirable  to  have  analytical  solutions  for 
simpler  geometries  to  facilitate  validation  against  benchmark  problems.8  The  traditional  way  to 
discretize  the  Schrodinger  equation  or  the  diffusion  equation  for  numerical  treatment  by  the  finite 
difference  method  is  to  employ  a  central  difference  formula  in  the  spatial  coordinates.  An  example 
of  this  is  the  implicit  Crank-Nicolson  scheme,  which,  for  the  diffusion  equation,  can  be  identified 
with  the  traditional  two- state  random  walk,  where  the  particle  is  endowed  with  two  directions  of 
motion  only.  Analytical  results  for  the  traditional  two-state  random  walk  with  perfectly  absorbing 
and/or  reflecting  boundaries  have  been  treated  in  number  of  works  including  Refs.  9  and  10. 
Separately,  the  case  of  the  telegraph  equation  with  partially  reflecting  boundaries  is  studied  in  Ref. 
11.  No  such  analytical  results  are  yet  available  for  the  4RW  model  and  it  is  the  purpose  of  the 
present  paper  to  provide  analytical  results  for  a  benchmark  initial-boundary-value  problem  in  1 
+ 1  space-time  dimension  for  the  model.  To  this  end,  we  extend  the  results  in  Ref.  7  by  considering 
two-point  boundary  conditions  for  the  4RW  model  and  derive  analytical  expressions  for  various 
transitional  probabilities.  Setting  aside  the  fact  that  the  4RW  model  has  a  physical  basis  in  the 
Feynman  chessboard  model  and  that  various  transitional  probabilities  are  related  to  the  solution  of 
the  continuous  Schrodinger  equation,  it  is  not  at  all  obvious  at  the  outset  from  the  governing 
difference  equations  that  an  analytical  solution  is  possible  for  the  said  boundary  value  problem, 
particularly  for  the  wavelike  solutions.  The  transform  approach  utilized  in  this  paper  will  reveal 
the  presence  of  the  discrete  Laplacian  operator  that  is  embedded  in  these  equations  and  will  clearly 
demonstrate  why  such  a  solution  is  still  possible,  while  paving  the  way  for  eigenfunction  expan¬ 
sion.  This  will  be  elaborated  in  Secs.  II  and  III.  The  most  general  case  of  partially  absorbing  and 
dissimilar  boundaries  is  considered  in  this  paper.  Results  for  the  special  cases  of  perfectly  absorb¬ 
ing  and  perfectly  reflecting  boundaries  are  also  provided  in  the  paper.  The  results  presented  here 
correspond  to  the  solution  of  the  discrete  form  of  the  diffusion  equation  as  well  as  to  the  real  and 
imaginary  parts  of  the  discrete  Schrodinger  equation. 

In  Sec.  II,  the  4RW  model  is  briefly  reviewed,  and  the  problem  under  investigation  is  defined. 
In  Sec.  Ill,  the  solution  to  the  4RW  model  subject  to  the  general  boundary  conditions  is  developed 
using  the  concept  of  generating  functions  and  eigenfunction  expansion.  Expressions  are  provided 
for  the  special  cases  of  perfectly  reflecting  and  perfectly  absorbing  boundaries  and  the  results  for 
the  latter  are  compared  to  numerical  simulations  obtained  by  directly  counting  paths.  Finally, 
conclusions  and  topics  of  further  research  are  given  in  Sec.  IV.  It  may  be  noted  parenthetically  that 
it  is  not  our  purpose  here  to  evaluate  other  various  models  that  have  been  put  forward  to  physically 
interpret  the  Schrodinger  equation,  a  topic  that  is  immensely  interesting  in  its  own  right. 


II.  4RW  MODEL 

For  a  particle  moving  on  a  discrete  lattice  and  subject  to  random  collisions,  the  transitional 
probabilities  considered  in  Ref.  7  at  the  discrete  space-time  point  (x=mA,t=se)  are  of  the  form 


(wi\  =  U(ex+e;1)  -(ex-etxx)\(wx 

\  w2)  2\  0  0  / \w2 
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E, 


<h 

<h 


J_/A‘ 

V2  W  Ex  )\q2r 


(2) 


where  q !  (m  A ,  s  6)  =  2  s/2\p  x  (m  A ,  s  e)  -p3(m\ ,  s  6)] ,  g2(m  A  ,se)  =  2s/2[p2(m\ ,  s  e)  -p4(m\ ,  s  6)  ] , 
Wi(mA,se)  =  [pi(mA,se)+p2(mA,se)+p3(mA,se)+p4(mA,se)],  w2(rnA,se)  =  [pi(mA,se) 

+p3(mA.,se)]-[p2(mA,se)+p4(mA,se)],  and  pIUL(mA,se)A,  jul=1,  ...  ,4,  is  the  probability  that  a 
particle  is  in  state  /z  at  the  space-time  point  (mA,se),  m=0,±l,±2,...,  s  =  0,l,....  The  particle 
changes  its  direction  of  motion  with  every  collision,  but  changes  its  parity  or  spin  at  every  other 
collision.  The  combination  of  two  directions  of  motion  and  two  states  of  spin  constitute  the  four 
states  in  the  model.  It  has  been  shown  in  Refs.  1  and  12  that  such  a  four- state  random  walk 
simultaneously  encompasses  the  diffusion  as  well  as  Schrodinger  equations.  The  quantities  nq  and 
w2  pertain  to  the  diffusion  process,  while  q2  and  ql  correspond  to  the  real  and  imaginary  parts  of 
the  Schrodinger  wave  function  in  the  discrete  case.  We  will  refer  to  (1)  as  the  diffusion  equation 
and  to  (2)  as  the  Schrodinger  equation  even  though  they  are  really  the  respective  discrete  coun¬ 
terparts  of  the  traditional  diffusion  and  Schrodinger  equations.  The  operators  Ex  and  Et  are, 
respectively,  the  spatial  and  temporal  advancing  operators  and  are  defined  mathematically  as 
E^1pJUi(mA,se)=pJUi[(m±l)A,se]  and  EtpfJi(mA,se)=p/il[mA,(s+  l)e].  It  is  assumed  in  Eqs.  (1) 
and  (2)  that  the  probability  that  a  particle  maintains  its  direction  at  the  next  time  step  remains  the 
same  as  the  probability  that  it  will  change  its  direction  at  the  next  time  step  and  that  the  probability 
of  remaining  at  the  same  location  at  the  next  time  step  is  zero.  If  the  number  of  right-going 
particles  is  the  same  as  those  going  to  the  left  at  time  t=0,  then  w2  =  0  and  Eq.  (1)  reduces  to  the 
simpler  equation 


E>,  =  Dxw  | ,  (3) 

where  DX=(EX+E~l)l 2  is  the  discrete  averaging  operator.  The  averaging  operator  in  (3)  owes  its 
existence  to  the  presence  of  the  Laplacian  operator  in  the  continuous  diffusion  equation  dw}/dt 
=Dld2wi/ dx2  when  the  spatial  and  temporal  step  sizes  are  subject  to  the  condition  A=  V2 eDx.  As 
such,  most  of  the  well-posed  issues  that  pertain  to  the  continuous  case13  will  be  carried  over  to  the 
discrete  case.  In  particular,  Eq.  (3)  will  be  well  posed  with  two-point  Robin  type  of  boundary 
conditions  and  the  solution  will  exist. 

The  difference  equations  (1)  and  (2)  are  assumed  to  be  valid  in  the  region  0 <m<€  and  s 
>0  and  they  are  supplemented  by  an  initial  condition  at  ^  =  0  and  boundary  conditions  at  m 
=  0,€.  We  will  adopt  the  abbreviation  v(m,s )  to  denote  the  discrete  function  v(mA,se).  The 
boundary  conditions  we  are  interested  in  are  of  the  form 

Pj(0,s)  -  axpj{\,s)  =  0  and  pj(€,s)  -  oi2pj{t  -  l,s)  =  0,  j=  (4) 

where  the  constants  a1  and  a2  are  assumed  to  be  real  and  positive.  These  are  the  discrete  versions 
of  the  general  Robin  type  of  boundary  conditions  for  the  continuous  case.  The  case  of  purely 
absorbing  boundaries  at  m  =  0,€  is  characterized  by  cq  =  0,  i- 1,2,  while  the  purely  reflecting  case 
is  characterized  by  a{  -  1 ,  i=  1,2. 9,14  The  general  case  corresponds  to  partially  absorbing  and 
partially  reflecting  boundaries  with  different  degrees  of  absorption  at  the  two  ends.  Our  interest  is 
to  obtain  analytical  solutions  to  (1)  and  (2)  on  a  discrete  space-time  lattice  (m\,se)  subject  to  the 
boundary  conditions  in  (4).  In  contrast  to  the  diffusion  equation  (3),  it  is  not  clear  at  the  outset 
whether  a  solution  will  exist  for  (2)  under  the  boundary  condition  (4),  setting  aside  the  fact  that  it 
is  tied  to  the  Schrodinger  equation.  This  is  because  of  the  presence  of  the  spatial  shift  operators 
that  are  neither  symmetric  (as  in  the  operator  Dx )  nor  asymmetric  (as  in  an  operator  of  the  form 
VX=[EX-E~1]).  Recall,  for  instance,  that  the  diffusion  equation  with  a  drift  term,  whose  discrete 
counterpart  will  have  neither  a  symmetric  nor  an  asymmetric  spatial  operator,  will  not  always  have 
a  solution  even  with  Neumann  type  of  boundary  conditions.  However,  we  will  demonstrate  in  Sec. 
Ill  that  the  temporally  transformed  equation  corresponding  to  (2)  will  indeed  contain  the  averaging 
operator  and  the  existence  question  will  be  set  to  rest.  Because  of  the  linearity  of  Eqs.  (1)  and  (2), 
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a  convenient  solution  can  be  obtained  by  using  generating  function  and  integral  transform  tech¬ 
niques  as  outlined  in  Refs.  7  and  15. 

III.  SOLUTION  BY  GENERATING  FUNCTIONS 

In  the  following,  we  assume  complete  symmetry  between  the  right-moving  and  left-moving 
particles  so  that  w2  =  0.  Consider  a  function  v(mA,se)  and  its  temporal  transform  d(mA,z)  as 
defined  in  Ref.  7, 


v(mA,z)  =  2  v{mA,se)zs  =  Tv.  (5) 

5=0 

The  quantity  d(mA,z)  may  be  thought  of  as  the  discrete  version  of  the  Laplace  transform  of 
v  and  is  simply  referred  to  as  the  z-transform.  The  inverse  relation  is  defined  as 


v(mA,se) 


v(mA,z) 


z 


5+1 


dz  =  Tlv, 


(6) 


where  Cz  is  a  closed  contour  around  the  origin  in  the  complex  z-plane  that  encloses  only  the 
singularities  at  the  origin,  /=  V— 1,  and  the  symbol  T  denotes  the  temporal  transform.  The  trans¬ 
formed  variable  v(mA,z)  is  also  referred  to  as  the  generating  function  within  the  random  walk 
community.9  Applying  the  temporal  transform  to  (2)  and  (3)  and  making  use  of  the  shift  property 
of  the  T transform  (7[£’?u1]=z_1[u1(mA,z)-i;1(m,0)])  and  carrying  out  some  simplifications,  we 
arrive  at  the  following  equations  for  various  generating  functions  in  terms  of  the  initial  conditions: 


[1  -  zDx]wi(mA,z)  =  W](m,0), 


(7) 


[1  -  ^2zDx  +  z2]qi(mA,z)  =  qi(m,0)  - 


z 

-j^Ex[q2(m,0)  +  <?i(m,0)], 


(8) 


[1  -  \f2zDx  +  z2]q2(mA,z)  =  q2(m,0)  -  -^Exl[q2(m,0)  -  <7i(m,0)].  (9) 

The  characteristic  operators  that  appear  on  the  left  hand  sides  of  (7)-(9)  are  generic  to  the  diffu¬ 
sion  and  the  Schrodinger  equations  under  the  4RW  model  and  are  seen  to  be  completely  indepen¬ 
dent  of  the  boundary  conditions.  A  remarkable  feature  of  the  spatial  dependence  of  these  operators, 
which  is  not  entirely  evident  in  the  initial  Eq.  (2)  describing  the  transitional  probabilities,  is  that 
they  all  involve  only  the  averaging  operator  Dx  that  is  the  discrete  counterpart  of  the  Laplacian 
operator  in  the  continuous  case.  Such  a  relation  has  already  been  alluded  to  in  Sec.  II.  Both  the 
diffusion  and  the  Schrodinger  equations  contain  the  Laplacian  operator  as  far  as  the  spatial  vari¬ 
ables  are  concerned,  and  the  transformation  from  the  continuous  case  to  the  discrete  case  for  a 
given  order  of  accuracy  is  not  unique.  Employing  a  central  difference  formula  for  the  spatial 
operator  will  lead  to  Crank-Nicolson  type  of  discrete  equations,8  which,  like  the  4RW  model, 
result  in  a  spatially  second  order  accurate  schemes.  The  important  point  to  note  from  (9)  is  that  the 
unknown  variable  contains  only  the  averaging  operator  Dx  and  that  any  other  asymmetries  that 
arise  from  Ex  or  E~l  alone  are  contained  only  on  the  right  hand  side,  operating  on  the  known  initial 
conditions.  A  formal  solution  to  Eqs.  (7)-(9)  can  be  affected  by  using  the  inverse  relation  (6)  and 
evaluating  the  integrals  in  the  complex  z-plane  after  expanding  the  reciprocal  of  the  characteristic 
operators  in  a  Maclaurin  series.  [Recall  that  the  contour  integral  in  the  inverse  operator  in  (6)  is  a 
closed  loop  of  vanishing  size  around  the  origin].  The  procedure  is  similar  to  that  outlined  in  Ref. 
7  and  will  involve  Chebyshev  polynomials  of  the  second  kind  for  the  wave  functions  qx  and  q2 . 
The  result  is 
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wx(m,s)  =  DsxWi(m,  0) , 


(10) 


D, 


A 


<h(m,s)  =  Us[  ~r=  )qi(m, 0)  -  -nUsA  ~j=  ^[?2(m,0)  +^1(m,0)], 


(11) 


^2(m,s)  =  Us 


D* 

A 


q2(m,0)  -  -\=U. 


Js- 1 


D, 

8 


E-x\q2(m,0)-qi(m,0)l 


(12) 


where  Us(x)  is  Chebyshev  polynomial  of  the  second  kind  of  order  s  and  argument  v.  The  formal 
solution  given  in  (10)-(12)  are  general  enough  and  are  valid  for  all  appropriate  boundary  condi¬ 
tions.  The  presence  of  the  operator  Dsx  in  the  diffusion  variable  wfm^s)  is  not  surprising  at  all. 
Indeed,  in  free  space,  (10)  directly  generates  the  characteristic  function  cos*  6  of  the  probability 
Wi  when  wfm,  0)  is  expanded  in  a  Fourier  series  with  transform  variable  0.  By  the  same  token, 
the  appearance  of  the  polynomials  Usf)  and  Us_f-)  is  intrinsic  to  the  Schrodinger  equation  in  the 
4RW  model,  as  already  indicated  previously.  Note  that  the  operators  Dx  (or  any  power  of  it)  and 
Ex  commute  and  the  order  of  these  terms  on  the  right  hand  sides  of  (11)  and  (12)  is  not  important. 
To  complete  the  solution,  we  must  now  expand  the  unknown  functions  in  terms  of  eigenfunctions 
of  the  Dx  operator  that  are  consistent  with  the  boundary  conditions  at  m  =  0,€.  In  free  space,  the 
appropriate  eigenfunctions  are  plane  waves  with  a  continuous  wavenumber  as  adopted  in  Ref.  7 
(or  said  in  other  words,  the  unknown  function  is  represented  in  terms  of  its  Fourier  transform  or 
series). 

For  the  boundary  conditions  indicated  as  in  (4),  we  seek  an  expansion  of  a  spatial  function 
v{m)  in  terms  of  exponential  and/or  trigonometric  functions.  An  exponential  function  of  the  form 
v0(m)  =  rm  is  a  valid  eigenfunction  provided  that  the  base  r=a\l  =  a2.  Clearly  this  is  only  possible 
in  the  special  case  of  ^ala2--=  ag=l,  where  ag  denotes  the  geometric  mean  of  a1  and  a2.  In  such 
a  case,  Dxu0(m)  =  0.5(a1  +  a2)v0(m)  :=  aa-v0(m),  where  aa,  being  the  arithmetic  mean  of  al  and 
a2 ,  is  the  eigenvalue  pertaining  to  u0(m).  When  a^l,  such  an  exponential  function  will  not  exist. 
We  will  denote  the  presence  of  this  exponential  function  by  employing  the  Kronecker  symbol  8la  . 
A  harmonic  function  of  the  form 


Uj(m )  =  Aj  sin (kjm)  +  Bj  cos (/ym),  7  =  1,2,...  (13) 

is  also  a  valid  eigenfunction  provided  that 

Bj(  1  -  a1  cos  kj )  =AJa1  sin  kj  (14) 

with  the  spatial  frequency  kj  given  by 

sin (kj£)  -  2 aa  sin [kj(£  -  1)]  +  a2g  sin [kj(£  -  2)]  =  0.  (15) 

Equations  (14)  and  (15)  are  obtained  by  enforcing  the  boundary  condition  (4)  at  the  two  ends.  A 
trivial  solution  of  Eq.  (15)  is  kj- 0.  There  will  be  a  total  of  (€-1)  nontrivial  solutions  of  (15),  thus 
constituting  a  total  of  €  distinct  eigenfunctions  [including  the  function  to  represent  the 

solution  of  (9)  and  (12)  at  the  points  m= 0,  ...,€.  Making  use  of  (14)  in  (13)  allows  us  to  extract 
a  bare  eigenfunction  Vj(m)  (i.e.,  without  the  coefficient  Aj  and  other  common  factors)  in  the  form 

Vj(m)  =  sin  (kjm)  -  ax  sin  [kj(m  -  1)],  (16) 

with  the  spatial  frequency  set  by  (15)  for  a  given  aq,  a2 ,  and  It  can  be  easily  verified  that 
DxVj(m)  =  cos  kj  •  Vj(m)  so  that  the  eigenvalue  of  Vj{m)  with  respect  to  the  averaging  operator  is 
cos  kj.  Furthermore,  the  eigenfunctions  Vj{m)  with  distinct  kj  as  well  as  Vj(m)  and  u0(m)  are 
mutually  orthogonal,  i.e., 
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€-1 

2  Vi{m)vj(m)  =  0,  ij  =  1,  ...,€- 1,  i  ^  j, 
m=l 


(17) 


€-1 

2  Vj(m)v0(m)  =  0,  j=  1, ...  1. 

m=  1 


(18) 


The  following  normalization  results  can  also  be  readily  established  for  the  functions  v0(m)  and 
vj(m): 


€-1 

X  vl(m)  = 

m=l 


-2(€-l)  _  i 


1-Off 


l/o0 


(19) 


and 


€-1 


X 

m=l 


sin  (k;£) 


sin  kj 


ax  cos[&;(€  -  2)]  -  -(1  +  acf)cos[kj(£  -  1)] 


€(1  +  a\) 

+ - -  cos  kj 


-  a\  sin 2[kj(£  -  1)]  :=  l/dj. 


(20) 


Gathering  all  of  the  above  results,  an  arbitrary  function  v(m)  satisfying  the  boundary  condition  (4) 
will  admit  a  spectral  representation  of  the  form 


€-1 

v(m)=A0v0(m)S1  +  X  AjVjim),  (21) 

7=1 

where  the  coefficients  A0  and  Aj  of  the  basis  functions  v0(m)  and  Vj(m)  can  be  obtained  in  terms 
of  the  function  v(m)  via  the  orthogonality  conditions  and  normalization  relations  (17)-(20), 


€-1 

A0  =  a0^  v0(m)v(m),  (22) 

m=  1 


€-1 

Aj  =  aj 2  Vj(m)v(m ),  j  =  1,  ...,€-  1.  (23) 

m=l 

With  this  spectral  representation,  it  is  clear  that 

€-1 

Dxv(m)  =  A0aa  •  u0(m)  ^  +  2  A;-  cos  kj  •  uy(m) .  (24) 

g  7=1 


Hence,  the  presence  of  Dx  in  the  spatial  domain  is  accounted  for  by  multiplying  the  spectral  basis 
function  term  by  its  respective  eigenvalue.  In  particular, 


€-1 

Dsxv(m)  =  A0asa  •  v0(m)  <5^+2  ^  cos*  kj  •  u7(m)  (25) 

*  7=1 


and 


t/^julm)  =  /l(,t/s(^)  •  L.°(m)^+  X  '  vAm)-  (26) 

Having  laid  out  the  formulation  for  the  general  case,  we  will  now  consider  two  interesting  special 
cases  which  will  merit  a  separate  discussion. 
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A.  Totally  absorbing  boundaries 

In  this  case  a1  =  a2= 0,  which  precludes  the  exponential  solution  v0(m).  Equation  (15)  implies 
a  solution  kj=irj/£,  j=  1 , ,€- 1.  The  eigenfunctions  Vj(m)  will  be  reduced  to  sin(7 rjm/€)  and 
cij=2/£  from  (20).  In  this  case,  one  gets  the  discrete  sine  transform16  representation  for  a  function 
satisfying  Dirichlet  boundary  conditions, 


€-1 

v(m)  =  2  Aj  sin  I 
7=1 


Aj=~^j  u(m)sinl 


'  m=  1 


jirm 

— 


7=1, 


(27) 


B.  Totally  reflecting  boundaries 

For  totally  reflecting  boundaries  at  both  ends,  a1  =  a2=  1,  which  dictate  that  ag=aa=l.  Equa¬ 
tion  (15)  indicates  a  solution  kj=jn/  (€-1) ,  y  =  1 1 .  The  exponential  function  now  reduces 
to  a  constant  function  v0(m)=  1.  The  normalization  constants  a0  and  aj  are  obtained  from  (19)  and 
(20)  as 


l/a0  =  ((  -  1),  l/aJ  =  2(€-l)sin2(^^y).  (28) 

An  unknown  function  v(m)  satisfying  perfectly  reflecting  conditions  at  the  two  ends  can  then  be 
expressed  as 


v{m)  =  A0  +  2  Ai  sin 
7=1 
€-1 

:=  A0  +  2  Cj  cos 


7 Tjm 


€  -  1 

- 1) 

€-1 


.  /  7rj(m-l) 
■  sin 


€  -  1 


=  A0  +  2  2 A  .-  sin 


7=1 


nj 


2(€-l) 


cos 


fr/(m  -  {) 


€  -  1 


(29) 


with  the  coefficients  given  by 


2  €_1  2  €_1 

TT^W.  ^=——2  cos 

1  ~  1  m=  1  C  “  1  m=l 


KJ 


-  D 


€-1 


u(m),  y  =  l,  ...,€-  1. 


(30) 


Equations  (29)  and  (30)  constitute  the  discrete-cosine-transform  (Ref.  16)  representation  of  a 
function  satisfying  Neumann  type  of  boundary  conditions  at  the  end  points. 


C.  Solution  for  general  case 

The  results  given  in  (25)  and  (26)  will  now  be  used  to  determine  the  Green’s  function13  of  the 
diffusion  equation  (10)  and  the  Schrodinger  equations  (11)  and  (12).  The  linearity  of  the  governing 
equations  with  respect  to  the  initial  conditions  enables  the  solution  to  arbitrary  initial  conditions  in 
terms  of  this  Green’s  function.  To  this  end  we  consider  initial  conditions  of  the  form  PjJjn,  0)A 
=  PfJL 0<5^°,  l<m0<€-l,  where  E^=1Pm0=  1 .  This  initial  condition  corresponds  to  the  case  where 
the  particles  are  released  from  the  location  ra0  in  a  state  /z  with  probability  P^0.  Letting  Tj 
=  (JP10-JP30)/A  and  T2  =  (P20-P40)/ A,  it  is  dear  from  the  definitions  that  w1(m,0)  =  ^°, 
^1(m,0)  =  rl^°,  and  q2(m,0)  =  T2S^°.  Using  the  spectral  representation  given  in  (21)-(23)  for  the 
functions  qi(m, 0),  and  q2(m,0),  we  arrive  at 

€-1 

Wi(m,0)  =  a0v0{m0)v0(m)Sla  +  2  ap  jim^v  j{m) ,  (31) 

*  h  i 
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€-1 


(him, 0)  =  fior i v0(m0)v0(m) d'a  +  a-v j(m0)v Am) 


J=  i 


€-1 


^r2(»i,0)  =  aor 2Vo(m,0)v0(m)  Sla  +  T22  av j(m0)v  Am) . 


(32) 


(33) 


j=i 


Substituting  (3 1)— (33)  into  (10)  and  (11)  and  recalling  relations  (25)  and  (26),  we  arrive  at  the 
solution  for  the  general  boundary  conditions  as 


€-1 


wx(m,s)  =  a0alv0(m0)v0(m)Sla  +  2  aj  coss(kj)vj(mQ)vj(m) 


(34) 


7=1 


and 


q1(m,s)  =  a0v0(m0)S1a 


€-1 

+  2  CLjVjirriQ ) 
7=1 


X 


,  cos  k;\  ,  N  r2  +  r,  /  cos  , 

x/2  / i  m  ~  ~^2rUs~i[~ir  \Vj{pi+l) 


(35) 


c/2(m,s)  =  a0v0(m0)Sla 


=  ]v()(m)  -  —  Sl  US_A %=  Mm-  1) 


r2c/„  ^ 


€-1 

+  2  ajVj(m0) 

j=  i 


X 


£ 


„  ,cosk:\  ,  r2-r,  cosk:,  , 

r,£/,|  — ^  k(m)  -  ^  Mm-  1) 


£  S-'V  £ 


(36) 


D.  Solution  for  totally  absorbing  case 

We  will  now  provide  some  numerical  results  for  the  special  case  of  totally  absorbing  bound¬ 
aries  at  the  two  ends.  In  this  case,  Eqs.  (34)-(36)  will  be  reduced  to 


.^-1 


2  v  .  /  77/ »?o \  .  /  JTjm\  TTj 

sm^— jcos 


(37) 


?i(m,j)  =  ^2sin|^pj 
/ 


57 

C0S\  €  /  |  .  ^  Trj in  'j 
^  - = —  ,cinl - 1 


£ 


'  \  €  / 


r2  +  r, 


77, 


5-1 


/  nj 

COS  -  11  / 

W  /  .  /  TTj(m  +  1) 

sin1 


V2 


(38) 
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FIG.  1.  Comparison  of  analytical  solution  for  w1(mA,s€)A  with  that  obtained  from  counting  paths  in  random  walk 
simulations  for  totally  absorbing  boundary  conditions. 


q2(m,s) 


r2-r, 


U, 


s- 1 


Trj(m-  1)\ 
€  / 


(39) 


Because  the  Green’s  functions  satisfy  reciprocity  conditions,  it  is  permissible  to  interchange  m  and 
m0  in  the  above  expressions  without  changing  the  field  values.  The  solution  for  a  continuous  space 
can  be  obtained  by  carrying  out  the  same  limiting  process  as  outlined  in  Ref.  7.  For  example,  Eq. 
(37)  reduces  to  the  expression  for  the  probability  density  function  of  the  diffusion  process  with  a 
diffusion  constant  D  with  two  totally  absorbing  points  at  v=0  and  x-L  and  an  impulsive  initial 
condition  at  v=v0, 


p<x-‘> = zf  izr)’  m 

a  result  that  is  found  in  many  texts  including  Ref.  9.  The  results  given  in  (34)-(39)  are  exact,  and 
as  such,  no  validation  is  necessary.  Nevertheless,  it  is  interesting  to  compare  the  numerical  data 
they  generate  with  those  generated  through  direct  random  walk  simulations,  as  the  latter  will  more 
likely  be  employed  for  more  complicated  time-dependent  boundaries.  In  such  cases,  the  analytical 
results  developed  here  will  serve  more  as  a  validation  check  for  the  numerical  results  generated 
through  random  walk  simulations.  The  analytical  results  (labeled  “Ana”  in  the  figures)  and  the 
results  obtained  by  counting  the  fractional  number  of  paths  (labeled  “Num”  in  the  figures)  in  the 
4RW  simulations  with  initial  conditions  set  at  F,10=0.5  =  F,2o  from  the  location  m0  =  4  are  shown  in 
Figs.  1  and  2  at  ^=16  and  for  €  =  32.  For  reference,  the  solution  for  Schrodinger  equation  in  free 
space  is  also  shown  in  Fig.  2.  For  the  parameters  chosen  in  Fig.  2,  the  solution  with  boundaries 
will  differ  from  the  free-space  case  only  for  0<m<  12,  as  is  clearly  seen  in  the  figure.  A  large 
number  of  realizations  (about  108)  was  needed  to  arrive  at  numerically  converging  results,  par¬ 
ticularly  for  the  Schrodinger  equation  whose  solution  involves  difference  in  probabilities.  Both 
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FIG.  2.  (Color  online)  Comparison  of  the  exact  solution  of  Schrodinger  equation  with  that  obtained  through  counting  paths 
in  random  walk  simulations. 


figures  reveal  that  the  numerical  simulated  results  closely  mimic  the  analytical  formulas  derived 
here.  The  numerical  results  also  underscore  the  importance  of  a  good  numerical  random  number 
generator,  particularly  needed  for  large  s,  as  the  sample  size  grows  exponentially  with  s.  The  use 
of  entwined  paths  described  in  Ref.  17  is  a  possible  means  of  getting  around  this  difficulty, 
although  details  have  only  been  demonstrated  there  for  the  relativistic  case  of  bounded  particle 
speed. 

E.  Appropriateness  of  the  4RW  model  to  wave  propagation  problems 

It  might  be  worthwhile  to  comment  a  little  on  the  relationship  between  the  4RW  model  and 
the  continuous  parabolic  equation  for  which  the  former  is  being  targeted  in  numerical  computa¬ 
tions.  For  the  parabolic  equation  encountered  in  a  number  of  wave  propagation  problems  such  as 
in  underwater  acoustics,  radio  wave  propagation,  and  optical  wave  propagation  in  fibers,  the 
governing  equation  for  the  field  variable  ^  in  a  homogeneous  medium  with  time-harmonic  exci¬ 
tation  is  of  the  form 


dijj  i  d^ijj 
dt  2k0  dx2  ’ 


(41) 


where  the  independent  variable  t  denotes  the  axial  spatial  coordinate  and  the  variable  v  denotes  the 
lateral  spatial  coordinate,  k0=27r/X  is  the  wavenumber  in  the  medium,  and  X  is  the  wavelength  of 
the  time-harmonic  excitation.  Equation  (41)  describes  two-dimensional  wave  propagation  in  the 
t-x- plane  subject  to  the  approximation  that  all  waves  travel  within  an  angle  of  9=  ±  15°  about  the 
axial  direction.3  It  can  be  shown  that  the  maximum  step  size  A  in  the  lateral  direction  is  restricted 
by  A^  X/(2  sin  0max),  where  #max  is  the  maximum  angle  of  propagation  with  respect  to  the  f-axis. 
As  6  is  restricted  to  a  value  less  than  15°  for  the  parabolic  approximation  to  be  valid,  we  may  take 
sin  9max  ~  tan  #max~  A/  e.  We  then  get  the  approximate  relation 


776  ^  k0A2 .  (42) 

It  may  be  remarked  parenthetically  that  the  inequality  (42)  translates  to  the  condition  that  the 
Schrodinger  equation  is  valid  for  describing  particle  motion  for  nonrelativistic  particle  speeds  with 
the  normalized  speed  v/c^ tan  #max~0.27  for  #max=15°,  where  c  is  the  speed  of  light  in  free 
space.  In  the  4RW  model,  the  relation  between  the  step  sizes  of  the  form  6=k0A2  is  implied  and 
this  is  consistent  with  the  inequality  (42)  If  A  =  X/V2,  then  this  relation  implies  that  6=7rX.  Both 
of  these  values  are  well  within  the  maximum  values  permitted  by  the  parabolic  equation  approxi¬ 
mation  and  it  is  believed  that  the  4RW  model  constitutes  a  very  appropriate  discretization  scheme 
for  numerically  handling  the  parabolic  equation. 
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IV.  CONCLUSIONS 

Analytical  expressions  have  been  provided  for  the  transitional  probabilities  of  a  4RW  model 
on  a  lattice  constrained  in  space  with  dissimilar  boundaries  and  subject  to  partially  reflecting 
boundary  conditions.  Special  cases  of  perfectly  absorbing  boundaries  and  perfectly  reflecting 
boundaries  have  been  treated.  Solution  for  the  transitional  probabilities  evaluated  at  time  s  to  the 
diffusion  process  is  shown  to  involve  sth  power  of  the  averaging  operator  Dx ,  whereas  those  of  the 
Schrodinger  equation  involve  Chebyshev  polynomials  of  the  second  kind  of  order  s  and  s- 1  with 
argument  Dx/  ^2.  These  are  the  general  characteristics  of  diffusion  and  Schrodinger  processes  in 
the  4RW  model  irrespective  of  the  boundary  conditions.  Different  boundary  conditions  will  dictate 
different  choices  of  the  eigenfunctions  in  which  the  initial  probabilities  at  s  =  0  are  expanded.  The 
eigenfunctions  in  the  most  general  case  will  involve  exponential  function  as  well  as  harmonic 
functions.  The  exponential  function  will  only  exist  when  the  parameters  present  in  the  boundary 
conditions  satisfy  certain  relationship.  While  the  results  presented  in  this  paper  should  have  a 
significance  of  their  own  for  random  walk  with  dissimilar  boundaries,  it  is  hoped  that  they  will 
also  serve  as  benchmark  cases  for  numerical  stochastic  methods  designed  to  solve  more  compli¬ 
cated  situations.  Although  it  had  not  been  the  major  focus  of  the  current  paper,  numerical  imple¬ 
mentation  of  the  random  walk  solution  necessitates  the  availability  of  an  effective  random  number 
generator  to  sufficiently  populate  all  portions  of  the  sample  space.  This  is  particularly  critical  for 
field  evaluated  at  large  times.  Extension  of  these  results  to  higher  dimensions  and  to  situations 
with  a  potential  field  is  a  topic  worthy  of  further  study  and  will  be  taken  up  in  the  future. 
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Abstract — The  angular  correlation  of  received  fields  in  multi- 
path  environments  is  studied.  The  focus  is  put  on  two-dimensional 
(2-D)  cases,  and  the  frequently  used  uncorrelated  scattering  as¬ 
sumption  is  tested  through  full-wave  Monte  Carlo  simulations.  The 
results  show  that  this  assumption  is  valid  for  the  discrete  finite 
spectra  of  the  received  waves  when  the  scattering  objects  in  the  en¬ 
vironments  are  distributed  in  complete  randomness,  either  when 
they  are  surrounding  the  transmit/receive  regions  or  in  clusters 
away  from  them.  The  correlation  among  wave  components  from 
different  angles  increases  only  when  the  randomness  of  the  scat- 
terer  deployment  is  reduced. 

Index  Terms — Angular  correlation,  discrete  finite  spectrum, 
random  media,  uncorrelated  scattering  assumption. 


I.  Introduction 

THE  angular  spectrum  of  certain  electromagnetic  (EM) 
field  in  the  presence  of  a  multipath  environment  is 
the  complex  amplitude  of  the  plane  wave  components  that 
consist  of  the  field  arriving  from  different  angles  at  a  fi¬ 
nite  sized  receive  volume.  Its  properties  are  of  special  use 
to  the  determination  of  the  multipath  richness  of  the  field, 
which  in  turn  is  crucial  to  the  performance  of  either  the 
simpler  spatial  diversity  schemes  [1]  or  the  more  advanced 
multiple-input-multiple-output  (MIMO)  systems  [2] .  The  mul¬ 
tipath  richness  depends  not  only  on  the  angular  spread  within 
which  waves  are  incident  onto  the  receive  volume,  but  also 
on  the  correlation  between  different  angular  wave  components 
[3].  Generally  speaking,  the  wider  the  angular  spread  and  the 
less  correlated  the  angular  components  are,  the  more  multipath 
richness  there  is  in  the  received  field. 

The  uncorrelated  scattering  assumption  is  an  angular  corre¬ 
lation  model  frequently  adopted  by  researchers  for  random  mul¬ 
tipath  propagation  studies  [4]— [8].  The  term  “uncorrelated  scat¬ 
tering”  here  is  consistent  with  the  one  adopted  by  a  recent  study 
of  Kennedy  et  al.  [8],  and  it  assumes  that  the  scattered  wave  re¬ 
ceived  from  different  angles  are  completely  uncorrelated.  Most 
of  the  time  it  is  further  assumed  that  the  angular  spectrum  is  also 
zero  mean,  which  will  be  shown  later  in  this  study  to  be  actually 
redundant.  This  uncorrelated  angular  scattering  assumption  is 
different  from  the  wide  sense  stationary  uncorrelated  scattering 
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(WSSUS)  assumption  proposed  by  Bello  [4]  on  a  first  look. 
However,  it  was  pointed  out  in  [3]  that  the  WSSUS  assump¬ 
tion  can  eventually  lead  to  the  uncorrelated  angular  scattering 
assumption,  which  is  the  subject  of  interest  here.  The  mathemat¬ 
ical  tractability  of  this  model  makes  it  very  popular  in  the  liter¬ 
ature  when  closed  form  conclusions  and  qualitative  interpreta¬ 
tions  are  sought  after.  Although  most  researchers  agree  that  this 
assumption  is  only  a  convenient  mathematical  idealization,  little 
work  has  been  done  to  determine  its  validity  and  how  well  it  rep¬ 
resents  a  true  random  environment.  Therefore,  the  objective  of 
this  work  is  to  study  the  autocorrelation  of  the  angular  spectra  of 
scattered  EM  waves  in  the  presence  of  discrete  random  media, 
and  to  investigate  the  validity  and  applicability  of  the  uncorre¬ 
lated  scattering  assumption. 

A  majority  part  of  this  work  is  carried  out  by  solving  scat¬ 
tering  problems  that  represent  multipath  environments  through 
surface  integral  equation  based  full  wave  techniques  [9],  [10]. 
The  operation  frequency  is  chosen  to  be  3  GHz.  For  the  2-D 
systems  under  consideration,  the  scattering  objects  are  mod¬ 
eled  by  lossy  dielectric  cylinders  of  different  shapes.  The  com¬ 
plex  dielectric  constant  is  taken  to  be  that  of  bricks  at  3  GHz, 
2.7  +  70.034  [11],  which  can  be  used  to  approximate  building 
structures  in  outdoor  environments.  Different  or  multiple  di¬ 
electric  constants  could  be  used,  and  it  will  only  influence  the 
complexity  of  the  numerical  calculations.  When  either  the  real 
permittivity  or  the  conductivity  of  the  material  varies,  the  cor¬ 
responding  wavelength  changes,  which  requires  different  dis¬ 
cretizations  on  different  object  surfaces.  However,  it  is  believed 
that  the  final  conclusions  of  the  study  will  be  independent  of  spe¬ 
cific  values.  As  will  be  shown  by  the  numerical  results  presented 
later,  the  validity  of  the  uncorrelated  scattering  assumption  is  de¬ 
termined  by  the  randomness  of  the  scattering  object  placement. 
For  more  efficient  solution  of  the  problems,  the  fast  multipole 
method  (FMM)  is  used  [12],  [13],  which  can  directly  results  in 
the  angular  spectra  of  the  fields  in  a  finite  receive  volume.  To  ob¬ 
tain  the  correlation  information  between  different  angular  wave 
components,  the  Monte  Carlo  simulation  technique  is  utilized, 
and  the  number  of  realizations  that  consist  the  sample  pool  is 
about  500  for  every  case  studied. 

The  rest  of  the  paper  is  structured  as  follows.  Section  II  in¬ 
troduces  the  basic  concepts  and  mathematical  essentials  of  the 
study.  Some  important  practical  considerations  concerning  the 
angular  spectra  are  pointed  out  in  Section  III.  The  angular  cor¬ 
relation  of  the  scattered  waves  is  then  studied  through  full  wave 
Monte  Carlo  simulations,  and  the  obtained  results  are  presented 
and  discussed  in  Section  IV.  Finally  Section  V  concludes  this 
study.  Throughout  the  work,  a  2-D  time-harmonic  propagation 
scenario  is  assumed.  The  e~XUJt  time  convention  is  adopted, 
where  u  is  the  radian  frequency  of  the  waves. 
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(a)  (b) 

Fig.  1.  Errors  in  {am}  due  to  ill-posedness  of  numerical  computation,  (a)  Configuration  of  a  scattering  system,  (b)  Computed  coefficients. 


II.  Theory 

A.  Angular  Spectrum 

Consider  a  finite  obstacle  free  receive  volume  Vr  in  two  di¬ 
mensions,  and  assume  that  it  can  be  enclosed  by  an  imaginary 
circle  of  radius  Rr ,  which  is  also  obstacle  free.  It  is  well  known 
that  the  scattered  field,  which  satisfies  the  2-D  homogeneous 
Helmholtz  equation,  takes  the  form 

oo 

Mp)=  E  amJm,(k0p)eim * 

m=  —  c o 
Mr 

«  u™Jm(k0p)eim*  (1) 

m=-MR+ 1 

in  the  circular  area  irrespective  of  the  sources  and  scattering  ob¬ 
jects  outside  [14].  In  (1),  Jm(-)  is  the  Bessel  function  of  the 
first  kind  and  of  order  ra,  p  =  (p,  $)  is  the  local  cylindrical 
coordinate  of  the  observation  point  with  respect  to  the  center 
of  the  circle,  ko  =  (27t/A)  is  the  wavenumber  in  the  back¬ 
ground  medium  (A  is  the  wavelength),  and  is  the  expansion 
coefficient  for  the  rath  cylindrical  harmonic.  For  a  given  accu¬ 
racy  requirement,  the  series  can  be  truncated  into  a  finite  length 
Mr  =  crkoRR ,  where  the  factor  a  is  a  constant  number  slightly 
greater  than  unity  and  controls  the  accuracy  of  the  truncation. 
It  is  indicated  in  [15]  that  a  sufficient  condition  to  ensure  the 
convergence  of  the  truncated  series  is  a  >  e/2.  Notice  must  be 
taken  that  for  the  series  to  converge,  it  does  not  require  the  co¬ 
efficients  Oim  to  approach  zero  when  |ra|  — v  oo.  As  a  matter  of 
fact,  it  is  the  exponential  decay  property  of  high  order  |  Jm(&op)  | 
that  ensures  the  convergence  [12],  [15],  [16]. 


By  using  the  integral  representation  of  the  Bessel  func¬ 
tions  Jm(z)  —  (1/27t)  Jq77  ^zcose errn(°-{^/2)) dO,  (1)  can  be 
rewritten  as 


where 


i>s(p) 


a(9)eikoP  cos(0~^  d8 
ae(ff)eikapcos^-s)dd 


(2) 


OO 

a(0)=  —  E  (3) 

m=  —  oo 
,  Mr 

ae(0)=^  E  «-eim(9_(7r/2))-  (4) 

rri  —  -  M  r-\-1 


act  as  the  complex  amplitudes  of  the  plane  wave  components 
of  the  received  scattered  field,  and  will  be  distinguished  from 
each  other  in  this  study  as  the  true  and  effective  finite  angular 
spectra  of  the  field,  respectively.  It  is  seen  that  both  spectra  are 
continuous  functions  of  propagation  angle  0 ,  and  all  the  plane 
wave  components  are  considered  to  be  propagating  waves.  Fur¬ 
thermore,  the  angular  spectra  are  in  the  form  of  Fourier  series 
whose  coefficients  are  27r},  and  it  is  straightforward 

to  obtain  the  cylindrical  harmonic  coefficients  {am}  from  the 
Fourier  coefficients  of  a(6)  by 


Analogous  to  representing  a  periodic  time  domain  signal  in 
its  Fourier  series,  the  different  coefficients  {am}  in  (3)  and 
(4)  correspond  to  different  frequencies  consisting  of  the  an¬ 
gular  spectra.  Clearly,  the  finite  spectrum  is  obtained  from  the 
true  one  by  neglecting  all  the  high  frequency  contents,  i.e., 
assuming  =  0  for  |ra|  >  Mr.  Since  the  actual  values  of 
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Fig.  2.  Finite  angular  spectra,  (a)  Continuous  spectra,  (b)  Discrete  spectra. 


am  in  this  range  may  not  be  negligible,  the  finite  spectrum  is 
not  necessarily  a  close  approximation  of  the  true  one.  However, 
because  Jm(p)  decays  exponentially  when  n  >  p,  the  two 
spectra  produce  practically  the  same  field  in  the  receive  volume. 
Without  causing  any  confusion,  the  term  frequency  will  be  used 
throughout  the  rest  of  this  work  referring  to  the  frequency 
components  of  angular  spectra. 

The  true  angular  spectrum  of  the  scattered  field  is  meaningful 
if  all  orders  of  {am}  are  available  a  priori.  On  the  other  hand, 
given  the  field  in  the  finite  receive  volume  only,  it  is  extremely 
difficult,  if  not  impossible,  to  obtain  its  true  angular  spectrum, 
because  the  high  frequency  components  produce  practically 
zero  field  in  the  volume  and  the  inverse  problem  becomes 
ill-posed.  Physically,  a  finite  aperture  only  provides  limited 
angular  resolution  [17],  [18].  In  order  to  appreciate  the  higher 
angular  frequencies,  the  size  of  the  aperture  has  to  be  increased. 
Due  to  this  reason,  what  really  matters  for  the  scattered  field  in 
the  receive  volume  is  the  finite  spectrum. 

B.  Uncorrelated  Scattering  Assumption 

In  a  random  environment,  the  angular  spectrum  of  the  field  in 
a  finite  receive  volume  is  a  random  process.  The  autocorrelation 
function  of  this  process  is  usually  called  the  angular  correlation 
of  the  field 

R(6,0,)  =  £{a(d)a*(6')}  (6) 

where  £{•}  represents  expectation  and  *  represents  complex 
conjugate.  The  uncorrelated  scattering  assumption  assumes  that 
the  scattered  waves  from  different  angles  are  completely  uncor¬ 
related,  which  can  be  represented  mathematically  as 

£  {a(6)a*  (O')}  =  £  {a(0)}  E  {a* (O')}  ,  0^0'  (7) 

for  nonzero  mean  a(0).  By  substituting  (3)  into  both  sides  of 
(7),  it  is  easy  to  arrive  at 

£{<**<}  =  £{<*»}£«}  (8) 

for  any  m  and  n.  Then  using  (8)  with  m  =  n,  one  can  find  out 
that  the  variance  of  is  zero  for  any  order  m ,  i.e., 

£  j|a%  -  £{am}\2 1  =  0,  Vm.  (9) 


This  suggests  that  the  sequence  of  “random  variables”  {am}  are 
actually  constants  with  probability  1,  and  hence  a(8)  becomes 
a  deterministic  function,  which  is  contradictory  to  the  starting 
point  that  angular  spectrum  is  a  random  process.  The  only  case 
of  resolving  this  contradiction  is  when  the  spectrum  is  a  zero 
mean  process,  for  which  the  uncorrelated  scattering  assumption 
takes  the  form 

£  {a(0)a*  (O')}  =  P(0)8(0  -  O')  (10) 

where  P(0)  is  the  angular  power  spectrum  proportional  to 
£  j|a(0)|2j,  and  £(•)  is  the  Dirac  delta  function.  One  can 
easily  check  (see  the  Appendix)  that  the  sufficient  and  neces¬ 
sary  condition  for  (10)  is 

p27T 

£{ama*}  =  /  P{e)e-^m~n">ede  (11) 

Jo 

which  is  a  function  of  only  the  order  difference  ( m  —  n )  re¬ 
lated  to  the  (m  —  n) th  Fourier  coefficient  of  the  angular  power 
spectrum.  Now  one  can  clearly  see  that  the  uncorrelated  scat¬ 
tering  assumption  for  the  true  angular  spectrum  (zero  mean)  of 
the  received  field  is  equivalent  to  a  stationary  coefficient  corre¬ 
lation  of  its  Fourier  series  expansion.  It  is  also  worth  mentioning 
here  that  the  Fourier  coefficients  of  the  angular  power  spectrum 
completely  determines  the  spatial  correlation  of  the  field  as  sug¬ 
gested  in  [6] . 

From  this  alternative  form  of  the  uncorrelated  scattering  as¬ 
sumption,  it  can  be  deduced  that  the  expansion  coefficient  <Ym 
may  not  become  negligibly  small  for  large  orders  almost  surely. 
Otherwise  setting  m  =  n  in  (11)  will  lead  to  E  ||am|2  j  = 

£  ||^oo  |2|  =  0  for  any  finite  m,  which  suggests  a  zero  scat¬ 
tered  field  in  the  receive  volume.  However,  it  is  reiterated  here 
that  even  though  the  high  frequency  components  of  the  true  an¬ 
gular  spectrum  may  have  significant  amplitudes  when  the  uncor¬ 
related  scattering  assumption  is  valid,  they  are  of  no  practical 
relevance  to  the  finite  receive  volume  because  of  the  resulting 
negligibly  small  field  caused  by  higher  order  Bessel  functions. 

III.  Practical  Considerations 

Given  the  field  in  a  finite-sized  receive  volume,  only  the 
coefficients  {<rm}  roughly  in  the  order  range  \m\  <  Mr  = 
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Fig.  3.  £{0  ™/'  *  }  for  systems  where  51  cylinders  are  around  the  receive  volume,  (a)  Real  part,  (b)  Imaginary  part. 


(ctyRr/X)  can  be  correctly  computed  with  high  confidence. 
Noise  and  numerical  errors  can  easily  blow  up  the  higher  orders 
during  the  inverse  process  due  to  their  ill-posedness.  One  way 
of  obtaining  is  to  solve  the  corresponding  scattering 

problem  by  the  means  of  surface  integral  equations  [9],  [10] 
and  apply  the  fast  multiple  method  (FMM)  [12],  [13],  provided 
that  all  the  scatterers  are  separated  from  the  receive  volume  by 
at  least  a  distance  that  is  comparable  to  the  size  of  the  volume 
(to  roughly  ensure  the  physics  of  propagating  plane  waves). 
Fig.  1  shows  a  system  with  a  single  circular  dielectric  scattering 
cylinder  and  the  computed  coefficients  {am}.  The  radius  of 
the  cylinder  is  6A  and  its  dielectric  constant  is  2.27  +  iO.034. 
The  transmitter  (TX)  is  a  unit-amplitude  electric  line  source, 
and  the  receive  volume  (RX)  is  a  circular  region  whose  radius 
is  taken  to  be  1.5A  and  3.0A  for  two  cases.  It  can  be  seen  from 
Fig.  1(b)  that  for  both  receive  volume  sizes  the  computed  higher 
order  coefficients  deviate  from  the  smooth  curve  formed  by  the 
lower  ones,  which  is  a  manifestation  of  the  errors  incurred.  The 
positions  of  the  deviation  are  marked  by  two  ellipses,  which 
are  around  \m\  =  12  and  \m\  =  23,  respectively.  These  results 
suggest  that  it  is  unlikely  to  obtain  the  true  angular  spectra  by 
manipulating  the  fields  in  a  finite  receive  region. 

In  principle,  the  true  spectra  can  be  acquired  by  measurement 
provided  that  the  equipment  is  ideal,  i.e.,  the  receive  antenna 
has  an  infinitely  narrow  beam.  However,  any  practical  antenna 
only  has  a  receive  pattern  with  finite  beam  width,  which  makes 
it  difficult  to  resolve  the  fast  oscillating  components  of  the  an¬ 
gular  spectra.  Therefore,  it  is  still  difficult  to  get  the  true  angular 
spectra,  and  the  left-hand  side  (LHS)  of  (10)  cannot  be  evaluated 
directly.  To  test  the  uncorrelated  scattering  assumption,  one  can 
only  rely  on  the  stationarity  of  (11)  in  the  correctly  computable 
low  frequency  range  as  a  partial  verification. 

Due  to  this  inability  to  determine  the  true  angular  spectra,  it 
is  more  reasonable  to  focus  on  the  finite  spectra  in  (4),  which 
produce  practically  the  same  field  in  the  receive  volume  as  do 
the  true  ones.  Furthermore,  since  a  finite  spectrum  excludes  all 
the  higher  frequencies,  it  is  also  possible  to  approximate  the 


Fig.  4.  1  —  |  Ce  (0m ,  Qn )  |  for  systems  where  5 1  cylinders  are  around  the  receive 
volume. 


corresponding  integration  in  (2)  as  a  summation  of  finite  terms 
following  the  trapezoidal  rule 

Nr-  1 

Mp)  «  X  ae{6n)eik°?m<*-e^A6  (12) 

n=0 

where  6n  =  (2 utt/Nr),  n  =  0, 1, ... ,  Nr  —  1  is  a  prede¬ 
fined  uniform  angular  grid,  A 6  =  (27 t/Nr)  is  the  grid  size, 
and  Nr  =  2 Mr  is  the  grid  number.  The  grid  number  here  is 
taken  to  be  the  same  as  the  number  of  observable  frequency 
components  in  the  spectrum,  and  the  accuracy  is  assured  as  sug¬ 
gested  by  Rokhlin  in  his  development  of  the  FMM  [12]  and 
the  bandwidth  consideration  mentioned  by  Chew  [16].  The  grid 
size  AO  can  also  be  interpreted  as  the  angular  resolution  of  the 
receive  volume,  and  (12)  suggests  that  this  resolution  is  fine 
enough  to  reconstruct  the  field  everywhere  inside  the  volume. 
If  one  further  defines  an  =  ae(0n)A0,  the  sequence  {an,  n  = 
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0, . . . ,  Nr  —  1}  constitutes  a  discrete  finite  spectrum  of  the  re¬ 
ceived  scattered  field. 

The  finite  spectrum  does  not  satisfy  the  uncorrelated  scat¬ 
tering  assumption  because  it  is  easy  to  check  that  the  correla¬ 
tion  between  its  Fourier  coefficients  of  various  orders  is  not  sta¬ 
tistically  stationary.  Even  so,  it  is  still  possible  that  two  of  its 
angular  components  that  are  separated  far  enough  are  uncorre¬ 
lated  or  only  weakly  correlated,  which  mathematically  means 

S  {ae(0)al(0')}  «  0,  for  ©  <  \6  -  0'\  <  2tt  -  ©  (13) 

with  0,0'  E  [0,  27t].  Furthermore,  due  to  the  limited  angular 
resolution  associated  with  the  finite  receive  volume,  it  is  un¬ 
necessary  to  seek  for  correlation  between  angles  separated  less 
than  the  resolution  size  AO.  Therefore,  it  would  be  interesting 
to  test  if  (13)  is  true  for  ©  =  AO.  In  the  numerical  results  pre¬ 
sented  in  Section  IV,  the  angular  correlation  will  only  be  evalu¬ 
ated  among  the  discrete  wave  angles  {On}.  Moreover,  for  better 
illustration  of  the  results,  the  normalized  autocorrelation  coeffi¬ 
cient  is  considered 

Ce{0m,8n)  = 

UK(U  -  £{ae{6m)}\  fc(ew)  -  £K(fln)}]}  (14) 

\f\ajjhn)  -  £{ffle(0™)}|2  -  £{a*(#n)}|2 

Another  issue  worth  mentioning  is  the  relation  between  the 
true  and  finite  angular  spectra.  Compared  to  the  true  spectrum, 
even  though  the  finite  one  ignores  all  the  high  frequency  compo¬ 
nents  am  for  \m\  >  Mr,  which  might  actually  have  significant 
amplitude,  it  is  still  possible  that  the  two  resemble  each  other 
quite  well  (see  the  example  that  follows).  When  the  discrete  fi¬ 
nite  spectra  are  considered,  the  errors  come  not  only  from  the 
frequency  domain  truncation,  but  also  possibly  from  angle  mis¬ 
matching,  which  is  also  demonstrated  by  the  following  example. 

When  a  single  plane  wave  is  incident  upon  a  circular  re¬ 
ceive  volume  centered  at  the  origin,  the  true  angular  spectrum 
is  clearly  nonzero  for  only  one  angle,  i.e.,  a(Q)  =  8(6  —  Oi), 
where  Oi  is  the  incident  angle.  The  corresponding  coefficients 
are  {am  =  which  are  of  unit  magnitude  up  to  infi¬ 

nite  order.  Suppose  the  radius  of  the  receive  volume  is  1.5A,  the 


finite  spectrum  includes  Nr  =  24  low  frequency  components. 
Fig.  2  shows  the  continuous  and  discrete  finite  angular  spectra 
ae(0)  and  an  with  incident  angle  Oi  =  315°  and  0%  =  320°.  The 
ripples  in  the  continuous  finite  spectra  of  Fig.  2(a)  are  obviously 
caused  by  ignoring  the  high  frequency  components.  However, 
the  strong  peaks  in  the  figure  correctly  identify  the  incident  an¬ 
gles.  For  the  discrete  finite  spectra  in  Fig.  2(b),  when  Oi  =  315°, 
it  is  correctly  indicated  that  there  is  only  one  angular  compo¬ 
nent.  However,  when  Oi  =  320°,  all  the  other  components  are 
also  present.  This  is  because  Oi  —  315°  is  one  of  the  pre-defined 
angles  {0n}  while  Oi  =  320°  is  not.  The  mismatching  between 
the  incident  angle  and  the  uniform  grid  requires  all  angular  com¬ 
ponents  to  appear  as  a  compensation. 

IV.  Simulation  Results 

This  section  studies  a  number  of  different  scenarios  and  the 
validity  of  the  uncorrelated  scattering  assumption  is  tested.  In 
the  first  system  studied,  a  line  source  transmitter  and  a  circular 
receive  volume  of  radius  1.5A  are  separated  by  100A.  Dielectric 
scattering  cylinders  of  different  shapes  are  randomly  distributed 
in  a  circular  ring  around  the  receive  volume  in  a  nearly  uniform 
way,  with  the  only  requirement  being  that  the  scatterers  do  not 
overlap  with  each  other.  The  inner  radius  of  the  ring  is  5A  and 
the  outer  radius  is  65 A.  The  shapes  of  the  cylinder  can  be  cir¬ 
cular,  elliptical,  or  square,  but  they  all  have  the  same  perimeter 
67rA.  Fig.  3  shows  the  grayscale  plots  of  the  real  and  imaginary 
parts  of  the  correlation  8{arna^l }  when  the  number  of  cylinders 
is  51  (corresponding  to  a  scatterer  area  density  around  11.5%). 
It  can  be  seen  that  £{ama;*  }  is  only  a  function  of  the  differ¬ 
ence  between  the  two  indices  m  and  n.  This  stationarity  sug¬ 
gests  that  it  is  likely  that  the  true  angular  spectra  satisfy  the  un¬ 
correlated  scattering  assumption.  To  accentuate  high  correlation 
with  dark  grayscale,  Fig.  4  shows  the  amplitude  of  the  simulated 
complementary  angular  correlation  of  the  discrete  finite  spectra, 
i.e.,  1  —  | Ge(0m,  0n) |.  For  the  complementary  correlation  coef¬ 
ficient,  values  close  to  zero  indicate  high  correlation.  The  waves 
from  two  distinct  discrete  angles  are  seen  to  be  almost  uncorre¬ 
lated,  suggesting  that  (13)  is  valid.  The  few  weakly  correlated 
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Fig.  7.  System  with  24  cylinders  distributed  in  a  cluster  away  from  transmit  and  receive  volumes,  (a)  |  £{ a  [  | .  (b)  1  —  |  Ce  (0m ,  Qn )  | . 


spots  around  the  diagonal  line  may  be  due  to  the  errors  from 
the  numerical  calculation  or  insufficient  number  of  realizations. 
Different  number  of  cylinders  are  also  used  to  see  the  effect  of 
scatterer  density.  The  results  are  similar  to  those  in  Figs.  3  and 
4,  suggesting  that  the  uncorrelated  scattering  assumption  is  not 
strongly  affected  by  the  scatterer  density  of  the  system. 

Next,  one  cylinder  of  each  shape  (3  total)  are  used  to  form  a 
group,  and  such  groups  of  cylinders  are  distributed  in  the  scat¬ 
tering  region.  The  relative  positions  of  the  cylinders  within  a 
group  are  fixed,  but  the  groups  can  still  be  randomly  located 
and  orientated.  Fig.  5  shows  the  results  with  4  groups  (12  cylin¬ 
ders),  and  similar  observations  to  before  can  be  made.  For  this 
case,  each  group  can  be  viewed  as  one  complicated  scattering 
object. 

Then,  a  system  similar  to  the  first  one  is  considered,  but  the 
randomness  of  the  cylinder  placement  is  decreased.  Fig.  6  shows 
the  results  with  36  cylinders  surrounding  the  receive  volume. 
However,  among  these  cylinders,  only  half  of  them  are  ran¬ 


domly  distributed,  with  the  other  half  taking  fixed  locations 
and  orientations.  It  is  observed  that  although  the  discrete  an¬ 
gular  spectrum  still  looks  somehow  uncorrelated,  the  correla¬ 
tion  E{arna^ }  is  no  longer  a  stationary  function,  and  hence  the 
uncorrelated  scattering  assumption  is  violated. 

In  another  type  of  system,  24  scattering  cylinders  are  placed 
in  a  circular  cluster  of  40A  radius  separated  from  both  the 
transmit  and  the  receive  volumes.  The  cluster  is  in  the  direction 
of  45°  to  the  transmitter  and  135°  to  the  receiver,  the  same 
geometry  as  Fig.  1(a)  with  the  single  cylinder  replaced  by  a 
cluster.  The  simulated  results  are  presented  in  Fig.  7.  It  can  be 
seen  that  the  coefficients  {am}  satisfy  the  stationarity  condition 
less  as  well  as  the  previous  cases.  The  finite  discrete  angular 
spectrum  appears  to  be  correlated  in  the  angular  range  0°-270° 
but  uncorrelated  from  270°  to  360°,  which  is  about  the  angular 
range  of  the  actual  scattered  waves  from  the  cluster.  Physically 
there  should  be  no  angular  components  whatsoever  coming 
from  the  range  0°-270°.  It  is  believed  that  the  seemingly 
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correlation  in  this  range  is  due  the  errors  introduced  by  trun¬ 
cating  the  series  in  (4)  and  angle  mismatching  discussed  in  the 
previous  section.  Hence,  it  is  still  possible  that  the  uncorrelated 
scattering  assumption  is  valid  for  the  true  spectrum,  and  (13)  is 
true  for  the  discrete  finite  spectrum. 

Next,  scatterers  are  distributed  in  a  circular  ring  around  the 
transmit  volume  instead  of  the  receive  volume.  The  size  of 
the  ring  is  the  same  as  before.  Fig.  8  shows  the  magnitudes 
of  and  the  complementary  correlation  coefficient 

1  —  \Ce(0m:0n)\  with  51  cylinders  as  scatterers.  Observations 
similar  to  the  case  of  Fig.  7  are  also  made  here,  where  the 
angular  range  of  uncorrelated  finite  spectrum  is  correctly  iden¬ 
tified  at  around  —  20°-20°.  The  same  arguments  as  the  previous 
case  can  be  made. 

In  the  last  system  studied,  24  scatterers  each  are  randomly 
distributed  in  circular  rings  around  both  the  transmit  volume 
and  the  receive  volume.  The  results  are  shown  in  Fig.  9,  and 


the  wave  components  from  different  angles  are  observed  to  be 
uncorrelated.  A  simple  comparison  with  Figs.  4  and  8  reveals 
that  the  wave  behavior  of  the  received  field  in  this  system  is 
dominated  by  the  influence  of  the  scatterers  immediately  around 
the  receive  volume,  which  is  a  very  reasonable  outcome. 

V.  Conclusion  and  Discussion 

By  studying  the  angular  components  of  the  scattered  field 
in  a  finite-sized  receive  volume  in  the  presence  of  various 
scattering  environments,  the  validity  of  the  uncorrelated  scat¬ 
tering  assumption  is  investigated  numerically  through  rigorous 
full  wave  electromagnetic  simulation.  The  results  about  the 
stationarity  of  the  spectrum  coefficients  indicate  that  it  is  likely 
that  the  assumption  is  generally  valid  for  the  true  spectra  in 
systems  where  scattering  objects  can  take  their  positions  and 
orientations  completely  randomly,  no  matter  they  are  around 
the  receive  volume,  the  transmitter,  or  in  clusters  separated 
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from  them.  In  the  meanwhile,  the  results  of  the  correlation  co¬ 
efficients  among  the  discrete  plane  wave  components  directly 
suggest  that  the  discrete  finite  spectra  are  also  likely  to  be  un¬ 
correlated  when  the  scattering  object  placement  is  completely 
random.  However,  as  soon  as  the  randomness  of  the  scatterer 
placement  is  decreased,  the  assumption  is  observed  to  break 
down  in  the  sense  that  the  correlation  of  the  spectrum  expansion 
coefficients  is  no  longer  stationary.  Intuitively,  there  is  a  direct 
cause-effect  relationship  between  the  randomness  of  the  scat¬ 
terer  deployment  and  the  validity  of  the  uncorrelated  scattering 
assumption.  In  outdoor  environments  where  a  large  number  of 
objects  can  move  around  freely,  this  assumption  might  be  a 
good  model  to  describe  the  multipath  wave  propagation. 

The  case  when  the  receiver  is  moving  with  constant  velocity 
is  of  special  interest  when  properties  such  as  space-time  cross 
correlation  and  space-frequency  cross  correlation  are  desired 
[19].  For  such  cases,  the  field  as  a  function  of  the  receiver  loca¬ 
tion,  which  in  turn  is  a  function  of  time  t,  can  be  obtained  from 
(2)  as  ?/;*(/.)  =  JQ27r  a(0)elkopo  cos(°-4>o) ^dt  cos(e-a) ^  where 

(pG,  Qo)  is  the  initial  location  of  the  receiver,  Ud  =  kov  is  the 
maximum  Doppler  frequency  shift,  and  a  is  the  direction  angle 
in  which  the  receiver  moves  at  velocity  v.  Similar  analysis  can 
be  carried  out  and  the  same  conclusion  about  the  angular  cor¬ 
relation  of  the  wave  spectrum  can  be  reached,  provided  that  the 
receiver  does  not  move  out  of  the  receive  volume.  Outside  the 
volume,  the  starting  (1)  may  not  be  valid. 

Furthermore,  even  though  the  study  presented  here  is  only 
for  the  received  field,  the  transmit  side  can  be  included  simply 
by  considering  reciprocity  of  the  systems,  and  the  uncorrelated 
scattering  assumption  can  be  easily  extended  to  the  double-an¬ 
gular  domain.  Specifically,  the  field  in  the  receive  volume  as 
the  response  to  some  source  in  the  transmit  volume  can  be  ex¬ 
pressed  as  a  combination  of  plane  waves  arriving  at  the  receive 
volume  in  different  angles  of  arrival  (AOA)  due  to  a  series  of 
plane  waves  leaving  the  transmit  volume  in  different  angles  of 
departure  (AOD).  When  the  scattering  objects  are  distributed  in 
complete  randomness,  the  coefficients  of  all  these  plane  wave 
components  are  uncorrelated  pairwise.  This  extended  assump¬ 
tion  has  been  successfully  applied  in  a  recent  work  of  the  au¬ 
thors  to  develop  a  statistical  propagation  model  between  finite 
transmit  and  receive  volumes  [20] . 

Appendix 

Alternative  Form  of  Uncorrelated 
Scattering  Assumption 

It  is  proved  here  that  the  sufficient  and  necessary  condition 
for  the  uncorrelated  scattering  assumption  in  (10)  is  (11).  The 
necessity  is  shown  first. 

Necessity:  By  substituting  (3)  in  the  LHS  of  (10)  and  using 
the  Fourier  expansion  of  the  Dirac  delta  function  on  the  right- 
hand  side  (RHS),  (10)  becomes 

o°  f  /  ,  \  2  oo 

E  U)  E 

n=  —  co  L  x  '  m=  —  c o 

o°  r  i 

xein(n/2)  e-in9'  =  J_p^ynd  9\  (15) 
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Matching  the  Fourier  coefficients  on  both  sides  and  some  simple 
manipulations  will  lead  to 


o°  r  i 

P{6)  =  V  e*(m-n)0  (16) 

Z — J  Z7T 

m=  —  c o  L  -1 

which  represents  the  power  angular  spectrum  in  a  Fourier  series. 
Using  the  standard  way  to  calculate  the  Fourier  coefficients,  one 
arrives  at  the  final  form  of  (11) 

/•2tt 

£{ama*n}  =  g^(m— n)(7r/2)  /  P(e)e-i{-m-n^d9.  (17) 


This  completes  the  proof  of  necessity. 

Sufficiency:  The  LHS  of  (10)  can  be  represented  as 


/  \  2  co  co 

£{a(8)a*(8')}  =  (-)  £ 


Xeim(9-(-K/2))  e-in(9' -(tt/2))  '  (jg) 


Substituting  (17)  into  (18),  one  obtains 


£{a(9)a*(9')}  =  /  P{8 1)  —  ^  e-im9leim9 


—  einf>1 e~ine'  J  dS 


P{61)6{61  -  8)S{8 i  -  0')d0\ 


=  P{8)8{8-8') 


which  is  exactly  (10).  This  completes  the  proof  of  sufficiency. 
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Abstract — We  discuss  the  use  of  the  parabolic  equation  (PE) 
along  with  the  alternate  direction  implicit  (ADI)  method  in  pre¬ 
dicting  the  loss  for  three  specialized  tunnel  cases:  curved  tunnels, 
branched  tunnels,  and  rough- walled  tunnels.  This  paper  builds  on 
previous  work  which  discusses  the  use  of  the  ADI-PE  in  modeling 
transmission  loss  in  smooth,  straight  tunnels.  For  each  specialized 
tunnel  case,  the  ADI-PE  formulation  is  presented  along  with 
necessary  boundary  conditions  and  tunnel  geometry  limitations. 
To  complete  the  study,  examples  are  presented  where  the  ADI-PE 
numerical  results  for  the  curved  and  rough-walled  tunnel  are 
compared  to  known  analytical  models  and  experimental  data,  and 
the  branched  tunnel  data  is  compared  to  the  numerical  solutions 
produced  by  HFSS. 

Index  Terms — Alternate  direction  implicit  (ADI),  parabolic 
equation,  radiowave  propagation,  tunnels. 


I.  INTRODUCTION 


THE  alternate  direction  implicit  (ADI)  method  coupled 
with  the  vector  parabolic  equation  (PE)  has  previously 
been  shown  to  model  radiowave  propagation  in  straight  tunnels 
with  smooth  walls  [1].  However,  due  to  the  rapid  growth  of 
telecommunication  systems,  different  tunnel  environments  also 
need  to  be  studied.  Subway  and  underground  road  tunnels 
typically  curve  or  branch  out  into  side  tunnels  and  have  walls 
which  are  not  smooth.  These  tunnel  geometries  are  not  always 
well  described  by  analytical  models  and  accurate  numerical 
solutions  become  important.  In  real  tunnels,  it  has  been  shown 
that,  over  a  long  distance,  high  order  modes  are  heavily  at¬ 
tenuated  and  low  order  modes  dominate  [3].  When  tunnels 
are  treated  as  imperfect  waveguides,  these  fields  represent 
waves  which  travel  near  the  axis  of  propagation.  The  PE  can 
accurately  model  low  order  modes  for  electrically  large  tunnels 
[1],  [2].  The  standard  PE  is  an  approximation  of  the  Helmholtz 
equation  that  assumes  fields  travel  within  ±15°  to  the  axis  of 
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propagation.  The  PE  lends  itself  to  numerical  discretization 
more  easily  than  the  Helmholtz  equation  but  does  not  account 
for  backscattered  fields.  Furthermore,  the  slow  varying  nature 
of  low  order  modes  implies  the  use  of  implicit  finite  difference 
(FDM)  techniques  where  large  discretizations  along  the  axis  of 
propagation  are  allowed. 

The  Crank-Nicolson  method  is  an  unconditionally  stable  im¬ 
plicit  FDM  that  has  been  traditionally  used  to  solve  for  the 
vector  PE.  However,  the  Crank-Nicolson  method  can  also  be¬ 
come  computationally  intensive  when  dealing  with  fine  meshes 
or  when  a  large  number  of  propagation  steps  is  required  [1] .  The 
alternate  direction  implicit  (ADI)  technique  is  a  modification  of 
the  Crank-Nicolson  method  that  reduces  computational  labor  by 
solving  for  the  fields  one  dimension  at  a  time  [1],  [6].  The  trun¬ 
cation  error  introduced  by  the  ADI  modification  is  of  the  same 
order  as  the  error  already  introduced  by  the  Crank-Nicolson 
method.  Previous  work  has  shown,  for  modest  discretizations, 
the  ADI  and  Crank-Nicolson  solutions  are  nearly  identical  [1]. 

In  this  paper,  we  use  the  vector  PE,  following  the  formula¬ 
tion  of  Popov  [2],  to  solve  for  fields  in  specialized  tunnel  envi¬ 
ronments.  For  each  case,  we  briefly  discuss  the  ADI  formula¬ 
tion  as  well  as  the  boundary  conditions  used  to  characterize  the 
tunnel  wall.  In  Section  II  we  discuss  the  curving  tunnel  and  com¬ 
pare  the  ADI-PE  results  to  known  analytical  approximations 
and  published  experimental  data  [7].  In  Section  III,  we  study 
the  branch  tunnel  and  compare  our  ADI-PE  numerical  results 
to  the  numerical  results  obtained  using  HFSS  [13]  for  a  smaller 
sample  problem.  Finally,  in  Section  IV,  we  formulate  a  model 
for  tunnels  with  surface  roughness  and  compare  our  ADI-PE  re¬ 
sults  with  known,  experimentally  verified,  analytical  solutions. 

II.  Tunnels  With  Smoothly  Curved  Axis 

A.  Curved  Tunnel  Propagation  Model 

Let  us  consider  a  tunnel  with  a  curved  axis.  The  geometry  of  a 
typical  curved  tunnel  with  a  rectangular  cross-section  is  shown 
in  Fig.  1,  where  s  is  the  curved  axis,  or  range,  and  p(s)  is  the 
range  dependant  radius  of  curvature.  The  vector  PE  was  formu¬ 
lated  by  Popov  and  was  shown  to  accurately  model  electromag¬ 
netic  propagation  in  curved  tunnels  [2] . 

The  vector  PE  for  a  tunnel  with  a  smooth  curve  in  the  hori¬ 
zontal  plane,  as  formulated  in  [2],  (with  e^ut  time  dependence, 
where  u  is  frequency  and  t  is  time)  is 


3W  d2W  d2W 

r,  2  "k  O  2 

os  oxz  oyz 
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Fig.  1.  A  typical  tunnel  with  a  curved  axis. 


where  kQ  is  the  free  space  wave  number,  and  W  is  a  vector 
function  that  is  directly  related  to  the  transverse  electric  field. 
The  relationship  of  W  to  the  transverse  electric  field  is  given  by 


W=(Em,Eyfejk°s.  (2) 


where  x  =  mAx,  y  =  lAy  and  s  =  sn  =  nAs.  The  differ¬ 
ence  quotients  are  8xWmj  =  Wm+1)i  -  2 Wmj  +  and 

8yWm,i  —  U7,,! , / 1  2 Wjn . i  T  lim./— i  and 


A 


n 

m 


B 


jmAxk0As\ 
2p(sn)  ) 
jmAxk0As ^ 

2p(*n)  J  * 


(7) 

(8) 


Note  that  when  p(s)  -a  oo,  A ^  and  B ^  become  unity  and  (5) 
and  (6)  reduces  to  the  ADI  formulation  for  the  straight  axis  PE 
shown  in  [1,  equations  (22)  and  (23)].  The  ADI  equations  in 
(5),  (6)  represents  a  marching  technique  where  the  transverse 
electric  field  is  solved  step  by  step  within  the  tunnel  domain  [2]. 
Starting  with  the  known  initial  field  at  the  n  =  0  plane,  the 
fields  of  each  successive  plane  is  solved  in  consecutive  order,  at 
propagation  steps  of  As,  until  the  field  at  the  desired  range  is 
solved. 


where  Ex  and  Ey  are  the  x  and  y  components  of  the  electric 
field,  respectively.  The  vector  PE  is  formulated  using  asymp¬ 
totic  analysis,  where  it  is  assumed  that  A  jb  <C  1  and  b/p(s )  1 
(where  A  is  the  wavelength),  which  means  it  is  only  valid  for 
high  frequency  propagation  and  in  tunnels  with  smooth  axis  of 
curvature. 

Along  the  tunnel  wall,  the  impedance  boundary  condition  is 
enforced  and  the  transverse  fields  become  coupled,  as  shown  by 
(3): 


jy  =  j_  f Ux  ny  \  (  zs  0  \  f  nx  ny 
ko  \^ly  Tlx  J  \  0  Zs  J  \Tly  Tlx 

where  nx  and  ny  are  the  x  and  y  components  of  the  unit 
normal  vector  at  the  boundary  and  Zs  is  the  normalized  surface 
impedance  [2] .  For  a  wall  with  relative  permittivity  er  and  con¬ 
ductivity  (j0  (in  S/m),  we  use  the  grazing  angle  approximation 
for  surface  impedance  [1],  [5] 


Erc 


where  erc  =  er  —  jar,  and  <jr  =  G0jujE0 ,  is  the  complex  permit¬ 
tivity  and  relative  conductivity,  respectively.  The  discretizations 
along  the  x,  y  and  s  axes  are  represented  by  Ax ,  Ay,  and  As, 
respectively. 

A  Peaceman-Rachford  [6]  ADI  formulation  of  (1)  can  be 
summarized  by 
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(5) 


(6) 


B.  Curved  Tunnel  Field  Approximations 

Approximate  analytical  solutions  describing  fields  in  curved 
tunnels  with  rectangular  cross  sections  and  constant  curvature 
radii  are  well  known  and  discussed  in  [3],  [8].  At  high  fre¬ 
quency,  kQa  1,  and  for  large  radii  of  curvature,  p  a,  the 
fields  are  best  described  in  terms  of  Airy  functions  [3],  [8].  As 
shown  in  (3),  the  nx  and  ny  components  of  the  normal  vector 
vanishes  alternatively  on  the  vertical  and  horizontal  walls  in 
tunnels  with  rectangular  cross  sections.  As  a  result,  the  com¬ 
ponents  of  the  vector  function  W  become  decoupled  and  can 
be  solved  independently.  The  horizontal  and  vertical  polariza¬ 
tions  can  now  be  solved  separately  by  enforcing  the  decoupled 
impedance  boundary  conditions  on  all  four  walls.  Following  the 
derivation  shown  in  [2],  the  vertical  component,  V,  of  the  vector 
function  is  shown  to  be  (keeping  the  same  notation  as  [2]): 


V  (x,  y,  s)  =  tt(x)E(y)eP22S  (9) 

Q(:r)  =  C\  Ai(t)  +  C<iBi(t )  (10) 

E(y)  =  D[Ai(x0)  -f  D2Bi(x0)  (11) 

where  Ai  and  Bi  are  the  Airy  functions  of  the  first  and  second 

1  /9 

kind,  P22  is  the  propagation  constant,  t  =  tQ-\-  (2 k2Q  / p)  x,  tQ 
and  xQ  are  constant  eigenvalues  to  be  found  from  the  boundary 
conditions,  given  by 


jkQ  dx 


=  0, 


x=l za 


At 


dV 


jkQZs  dy 


y=0,b 


=  0. 
(12) 


For  large  curvature  radii,  the  vertically  polarized  fields  are 


n(x) 

Z(y) 


sin  ■ 


mry 


1,2,3... 


1,2,3... 


(13) 

(14) 


where  the  integers  to  and  n  physically  represent  field  variations 
along  the  x  and  y  axes  and  specifies  a  possible  mode.  The  mode 
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TABLE  I  TABLE  II 

Analytical  and  Numerical  MAFs  for  the  Curved  Tunnel  With  Analytical  and  Numerical  MAFs  for  the  Curved  Tunnel  With 

Rectangular  Cross  Section  for  950  MHz  Rectangular  Cross  Section  for  1.8  GHz 


f  =  950  MHz 

P  (m) 

Analytical  (dB/km) 

ADI  (dB/km) 

oo 

16.2 

16.0 

2533 

17.3 

18.1 

1267 

18.8 

18.5 

844 

20.3 

19.8 

633 

21.7 

20.8 

507 

23.2 

22.0 

f  =  1.8  GHz 

P  (m) 

Analytical  (dB/km) 

ADI  (dB/km) 

oo 

4.5 

4.5 

4800 

5.2 

5.2 

2400 

6.0 

6.2 

1600 

6.7 

6.7 

1200 

7.5 

7.5 

960 

8.2 

8.2 

attenuation  per  unit  length  can  be  found  from  the  complex  ex¬ 
ponent,  P22,  to  be  [2] 


az 


dB 

km 


:  4343A 2 Re  (  —  )  +  4343A 2 Re  (Zs) 
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A2P 


_(2a)J 


+ 


8(2a)3  (15— m27r2) 
37r2(mA)4p2 


t+  >  l 

tt «  -i 


(15) 


where 


and  q  =  (2 k^/(p  a))1^3  a.  The  £+  >  1  part  of  (15)  de¬ 


scribes  small  radii  of  curvature  and  the  <C  —  1  part  de¬ 
scribes  large  radii  of  curvature  (or  almost  straight  tunnels).  Mah¬ 
moud  and  Wait  [9]  numerically  calculated  the  attenuation  fac¬ 
tors  for  curved  waveguides  using  more  precise  transcendental 
equations.  As  shown  in  [2,  Fig.  3],  the  aysmptotic  equation  (15) 
is  in  very  good  agreement  with  Mahmoud  and  Wait’ s  numerical 
solutions.  We  will  use  the  aysmptotic  equation  to  validate  the 
numerical  simulations  of  the  ADI-PE. 


C.  Comparison  of  the  ADI-PE  to  Analytical  Solutions 

In  this  sub- section,  we  validate  the  ADI-PE  with  numerical 
examples.  We  consider  a  rectangular  tunnel  with  dimensions 
of  8  m  x  4  m  and  walls  with  relative  dielectric  constant,  er  = 
5.5,  and  conductivity,  cr0  =  0.01  S/m.  We  used  the  dominant 
m  =  1,  n  =  1,  mode  of  the  vertically  polarized  field,  as  shown 
by  (13)  and  (14),  as  our  initial  field.  The  discretizations  used 
for  the  950  MHz  case  were  As  =  3.73  A,  Ax  =  0.84  A  and 
Ay  =  0.70  A;  and  the  discretizations  used  for  the  1.8  GHz 
simulations  were  As  =  3.75  A,  Ax  =  0.80  A  and  Ay  =  0.71  A. 

The  mode  attenuations  per  unit  length,  or  mode  attenuation 
factors  (MAFs),  from  (15)  and  the  ADI-PE  simulations  are  tab¬ 
ulated  in  Tables  I  and  II  for  different  curvature  radii. 

As  shown  in  Tables  I  and  II,  the  numerical  MAFs  closely 
follow  the  aysmptotic  solutions  over  the  range  of  curvature  radii. 
Furthermore,  as  the  curvature  radii  decreases,  the  number  of 
reflections  from  the  walls  along  the  curved  path  increases  and 
so  does  the  loss.  The  percent  error,  defined  by 


%Attenuation  Error  = 


\MAFADI-ax(£)\ 


(17) 


Fig.  2.  Percent  error  of  the  MAFs  of  the  curved  waveguide  for  950  MHz  (dot, 
solid)  and  1.8  GHz  (circle,  solid)  for  As  =  3.75  A.  The  percent  error  of 
the  MAFs  for  950  MHz  (dot, dashed)  and  1.8  GHz  (circle,  dashed)  for  As  = 
1.875  A.  The  vertical  lines  show  the  region  where  t+  =  ±1  for  950  MHz 
(solid)  and  1.8  GHz  (dashed). 


is  shown  in  Fig.  2  as  a  function  of  D2 /Xp,  (where  D  is  the  length 
of  the  diagonal  of  the  rectangular  cross  section)  for  950  MHz 
(dot,  solid)  and  1.8  GHz  (circle,  solid).  The  figure  also  shows  the 
results  for  smaller  discretizations  (As  =  1.875  A)  with  dashed 
curves  to  show  that  the  solutions  are  convergent.  The  vertical 
lines  in  Fig.  2  show  the  region  where  t+  =  ±1  for  the  950  MHz 
(solid)  and  1.8  GHz  (dashed)  cases,  respectively.  The  parameter, 
D2/Xp ,  was  selected  from  [2]  because  it  is  a  wavelength  depen¬ 
dent  term  that  must  be  much  less  than  unity  to  ensure  accurate 
results.  Furthermore,  D2 /Xp,  is  the  Fresnel  number  and  is  one 
of  the  parameters  which  govern  the  diffraction  processes  in  the 
waveguide  [2].  As  Fig.  2  shows,  there  is  less  than  5%  error  in 
MAFs  over  the  entire  range  of  p  considered;  even  when  there 
are  sharp  bends  and  the  condition,  D2 /Xp  <C  1,  is  not  satisfied. 
Finally,  one  of  the  important  characteristics  of  the  curved  tunnel 
is  the  accumulation  of  the  field  near  the  concave  wall  at  y  —  —a 
(the  whispering  gallery  mode).  Fig.  3(a)  shows  the  dominant 
| Eli  |  field  as  defined  by  (13)  and  (14)  across  the  initial  s  —  0 
plane  at  1.8  GHz.  The  whispering  gallery  mode  feature  can  be 
seen  in  Fig.  3(b),  where  the  field  generated  by  the  ADI-PE,  at 
a  range  of  2000  m  for  p  =  2000  m,  has  an  accumulation  along 
the  x  =  — 4  m  wall. 

D.  Comparison  of  the  ADI-PE  to  Ray  Tracing  Simulations 

In  this  section,  we  compare  the  ADI-PE  simulations  to  the  ray 
tracing  simulations  shown  in  Wang  et  al.  for  a  curved  waveguide 
[7] .  The  geometry  of  the  waveguide  (top  view)  is  shown  in  Fig.  4 
and  is  comprised  of  two  straight  sections  and  a  curved  section. 
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x(m)  x(m) 

Fig.  3.  The  \  Vn  \  mode  at  (a)  s  =  0  m  and  at  (b)  5  =  2000  m  for  p  —  2000  m 
at  1.8  GHz. 


Curved  Section 


Fig.  4.  Geometry  of  the  curved  tunnel  with  straight  sections. 


The  axial  length  of  the  waveguide  is  400  m  and  the  length  of 
the  curved  section  is  200  m  with  a  radius  of  curvature  of  300  m. 
The  cross-sectional  dimensions  of  the  waveguide  are  8  m  x  6 
m.  In  [7],  the  transmitter  and  receivers  are  vertically  polarized 
half-wave  dipoles  operating  at  a  frequency  of  1  GHz.  The  trans¬ 
mitter  and  receivers  are  located  at  the  center  of  the  waveguide 
at  heights  of  3  and  1.5  m. 

To  establish  a  basis  of  comparison,  we  first  look  at  the  straight 
waveguide.  Fig.  5  shows  the  normalized  received  power  of  the 
ray  tracing  simulations  (solid)  and  the  ADI-PE  simulations 
(dashed).  The  ADI-PE  is  simulated  using  the  far  field  expres¬ 
sions  of  a  half-wave  dipole  in  free  space  placed  30  m  outside 
the  entrance  as  the  initial  field.  The  field  is  tapered  by  a  unit 
Gaussian  with  standard  deviation,  a  =  2.5  A  to  minimize  error 
associated  with  using  incorrect  field  values  at  the  walls  [1].  The 
dielectric  constant  and  conductivity  of  the  waveguide  walls  are 
taken  to  be  the  typical  values  of  5  and  0.01  S/m,  respectively. 
Simulations  show  that  over  the  range  of  acceptable  values  of 
dielectric  constants  (er  =  5  —  12),  there  is  little  change  in  the 
received  power.  The  simulations  are  done  with  discretizations 
of  Ax  =  0.67  A,  Ay  =  0.74  A  and  As  =  1.33  A.  As  the  figure 
shows,  the  ADI-PE  closely  models  the  nulls  in  the  received 
power  generated  from  ray  tracing  simulations.  The  ADI-PE 
results  are  vertically  offset  so  that  the  least  squares  fit  line  of  the 
ADI-PE  data  and  the  experimental  data  (in  the  curved  region) 
intersect  at  the  start  of  the  curved  region  (at  100  m). 

As  discussed  in  [1]  and  [3],  propagation  in  straight  tunnels 
can  be  characterized  by  a  near  and  far  zone.  In  the  near  zone, 
rays  propagating  at  large  angles  make  significant  contributions 
to  the  field  and  take  the  form  of  rapid  oscillations.  In  the  far 
zone,  these  rays  are  severely  attenuated  and  paraxial  rays  are 
dominant.  In  the  near  zone,  the  PE  is  not  accurate,  but  so  long 
as  the  low  order  modes  are  illuminated,  the  PE  will  be  accurate 
in  the  far  zone.  The  start  of  the  far  zone  is  determined  by  the 
size,  shape  and  frequency  of  the  tunnel  [3].  The  far  zone  can 
be  found  by  calculating  the  attenuation  constant  of  each  mode. 
However,  the  attenuation  constant  and  amplification  term  for 
each  mode  may  not  always  be  found  analytically.  If  we  summed 
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Fig.  5.  Normalized  received  power  from  ray  tracing  (solid)  and  ADI-PE 
(dashed)  along  the  axis  of  propagation  for  the  straight  tunnel. 
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Fig.  6.  Normalized  received  power  from  ray  tracing  (solid)  and  ADI-PE 
(dashed)  along  the  axis  of  propagation  for  the  curved  tunnel. 


the  equally  weighted  first  100  modes  for  the  rectangular  waveg¬ 
uides  using  the  propagation  constant  expressions  defined  by  [3, 
Eqs.  (54)  and  (55)],  we  can  see  that  the  region  before  500  m 
is  in  the  near  zone.  In  the  range  considered  here,  the  contribu¬ 
tions  of  the  high  order  modes  account  for  the  discrepancies  in 
the  ADI-PE  and  ray  tracing  results.  These  discrepancies  are  in 
the  form  of  rapid  oscillations  in  the  experimental  data  along  the 
axis  of  propagation.  Fig.  6  shows  the  normalized  received  power 
for  the  curved  waveguide.  As  shown  in  the  figure,  the  ADI-PE 
closely  tracks  the  received  power  of  the  ray  tracing  simulations 
within  the  curved  section. 

E.  Comparison  of  the  ADI-PE  to  Experimental  Data 

In  this  section,  we  compare  the  ADI-PE  to  experimental  data 
for  the  Lin-sen  subway  tunnel,  shown  in  [7,  Fig.  6].  The  Lin-sen 
subway  is  comprised  of  two  curved  sections  with  radii  of  cur¬ 
vature  of  455.68  m  and  354.74  m,  respectively,  separated  by  a 
straight  section.  The  straight  entrance  and  exit  sections  are  not 
considered  here  because  they  do  not  lie  within  the  same  hori¬ 
zontal  plane  as  the  rest  of  the  tunnel.  The  tunnel  cross-section  is 
approximately  rectangular  with  dimensions  of  6  m  x  8  m.  The 
transmitter  is  located  outside  the  tunnel  and  is  a  vertically  po¬ 
larized  Yagi-Uda  antenna  operating  at  942  MHz.  The  receiver 
is  placed  off  center  at  a  height  of  1.85  m  above  the  ground. 

Fig.  7  shows  the  received  power  from  measurements  (solid) 
and  the  ADI-PE  simulations.  As  before,  the  dielectric  constant 
and  conductivity  are  assumed  to  be  the  typical  values  of  5  and 
0.01  S/m.  In  this  case,  the  field  entering  the  horizontal  section  of 
the  subway  is  unknown  and  represents  a  possible  source  of  error. 
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Fig.  7.  Experimental  data  (solid)  and  ADI  results  with  the  HE2i 
(long-dashed),  HEX1  (short-dashed)  mode  and  the  half  wave  dipole  (dot-solid) 
as  the  initial  field  for  the  Lin-sen  subway  tunnel. 


However,  we  may  choose  a  low  order  mode  as  our  initial  field 
because,  at  large  distances,  only  low  order  modes  dominate  [1]. 
Simulations  done  for  various  initial  fields  show  that  the  m  = 
2,n  =  1  mode  of  the  field  (long-dashed),  as  defined  by  (13) 
and  (14),  produces  solutions  that  best  fits  the  experimental  data. 
The  discretizations  are  Ax  =  Ay  =  0.72  A,  and  A z  =  2.11  A. 
As  Fig.  7  shows  for  this  initial  field,  even  with  our  limitations, 
there  is  still  good  agreement  between  the  ADI-PE  solutions  and 
the  measured  results.  The  ADI-PE  models  the  nulls  and  overall 
trend  of  the  received  power.  The  figure  also  shows  the  field  in¬ 
tensity  when  the  m  =  l,n  =  1  mode  (short-dashed)  and  the 
far  field  expressions  for  the  half-wave  dipole  (dot-line)  are  used 
as  initial  field.  What  this  shows  is  that  the  field  intensity  is  very 
sensitive  to  the  initial  field.  This  is  because  we  are  dealing  with 
a  short  curved  section  (approximately  50  m),  and  high  order 
modes  make  significant  contributions  in  to  the  straight  section. 
As  in  the  previous  example,  we  are  not  operating  in  the  far  zone 
and  high  order  modes  make  significant  contributions.  The  high 
order  modes  are  represented  by  the  rapid  fluctuations  in  the  mea¬ 
sured  data.  The  PE  approximation  does  not  accurately  model 
these  modes  and  it  is  a  source  of  error.  These  results  were  ob¬ 
tained  in  a  matter  of  minutes  using  a  typical  PC  (1  GB  RAM). 
By  comparison,  a  typical  ray  tracing  code  would  require  another 
simulator  to  generate  the  initial  field  and  must  track  each  reflec¬ 
tion  and  diffraction  for  each  ray  and  can  become  computation¬ 
ally  intensive. 

III.  Branch  Tunnels 
A.  Branch  Tunnel  Model 

Let  us  now  consider  the  case  of  a  straight  tunnel  that  branches 
into  a  side  tunnel.  A  typical  branch  tunnel  geometry  is  shown 
in  Fig.  8(a)  and  (b),  where  the  main  tunnel  axis  is  shown  as  a 
solid  bold  line  and  the  branch  tunnel  axis  is  shown  as  a  long 
dashed  line.  The  branch  angle,  Ob,  is  the  angle  between  the  axes 
of  the  straight  and  branch  tunnels.  Fig.  8(a)  and  (b)  show  the  in¬ 
cident  and  reflected  rays  as  it  enters  the  branch.  The  short  dashed 
line  marks  the  input  plane  of  the  branch  tunnel.  The  grazing 


(a)  (b) 

Fig.  8.  The  incident  and  reflected  rays  entering  the  branch  tunnel  when  (a) 
'E  <  0b  and  (b)  xlf  >  6b. 


angle,  and  the  angle  of  the  ray  entering  the  branch,  a,  are 
also  shown.  In  order  for  the  PE  approximation  to  be  valid,  the 
branch  angle  must  be  small  enough.  More  precisely,  the  branch 
angle  must  be  less  than  30°  to  ensure  reflected  rays  entering  the 
branch  are  within  our  PE  limit  of  ±15°.  Also,  considering  the 
rays  diffracted  from  the  comers  of  the  junction,  we  can  arrive 
at  a  much  more  stringent  requirement  of  Ob  <  15°.  Diffracted 
rays  entering  the  main  tunnel  at  angles  greater  than  15°  will  be 
weak  when  compared  to  reflected  rays  and  will  experience  se¬ 
vere  attenuation  after  the  tunnel  junction. 

As  in  the  previous  section,  we  solve  the  vector  PE  shown 
in  (1)  with  p(s)  — ►  oo.  The  slope  of  the  branching  wall  is 
modeled  using  a  staircase  approximation  (see  Fig.  9(a))  and  the 
impedance  boundary  condition  is  enforced  on  all  four  walls  as 
outlined  in  [1]  and  [2].  The  fields  along  the  planes  marking  the 
entrance  of  the  main  tunnel  (line  C)  and  the  branch  tunnel  (line 
B)  are  solved  simultaneously  and  then  used  as  the  initial  fields 
for  the  two  separate  diverging  tunnels. 

B.  Comparison  of  ADI-PE  Results  to  HFSS  Simulations 

In  this  section  we  validate  the  ADI-PE  branch  model  with  a 
numerical  example.  We  simulate  a  0.9  m  x  1.0  m  rectangular 
tunnel  with  a  branch  angle  of  15°  and  operating  at  a  frequency 
of  900  MHz.  The  initial  field  is  a  unit  strength  Gaussian  field 
source  in  the  far  zone.  The  source,  with  standard  deviation  of 
0.75  A,  is  centered  5  m  before  the  tunnel  junction  (only  the  re¬ 
gion  near  the  tunnel  junction  is  shown  in  Fig.  9).  This  means  we 
are  only  using  the  lowest  order  modes  as  our  initial  field  [1],  [3]. 
The  fundamental  mode  propagates  near  our  PE  limit  at  an  angle 
of  14°  with  respect  to  the  axis  of  the  main  branch. 

The  ADI-PE  simulations  are  done  with  discretizations 
(within  the  tunnel  junction)  of  Ax  =  0.15  A,  Ay  =  0.14  A  and 
Az  =  0.13  A.  The  cross-sectional  coordinates  are  indicated 
by  x  and  y,  while  the  axial  coordinate  in  the  main  tunnel  is 
denoted  by  the  z-axis.  The  z-axis  discretizations  are  made 
small  within  the  tunnel  junction  to  ensure  small  step  sizes  for 
the  staircase  approximation.  Outside  the  junction  region,  dis¬ 
cretizations  along  the  axis  of  propagation  can  be  made  as  large 
as  a  few  wavelengths  [1].  To  validate  our  results,  we  compare 
our  solutions  with  HFSS  [13]  and  plot  the  magnitude  of  the  Ey 
field  along  the  main  and  branch  tunnel  axes  in  Fig.  9(b).  The 
HFSS  simulations  use  radiation  boundary  conditions  to  termi¬ 
nate  the  tunnel  and  symmetry  planes  to  reduce  computational 
labor.  HFSS  is  a  full  wave  simulator  and,  unlike  the  ADI-PE, 
solves  for  backscattered  waves  as  well  as  waves  traveling  in  the 
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Fig.  9.  (a)  Geometry  of  the  branch  tunnel,  (b)  The  axial  field  intensity  of  the  main  tunnel  from  ADI-PE  (solid),  HFSS  (dashed),  the  branch  tunnel  from  ADI-PE 
(asterisk)  and  HFSS  (dashed). 


forward  direction.  The  backscattering  is  seen  as  fluctuations 
in  the  axial  field  in  Fig.  9(b)  near  the  diverging  tunnels.  The 
dielectric  constant  and  conductivity  of  the  tunnel  walls  is  5  and 
0.1  S/m,  respectively.  A  high  conductivity  is  chosen  so  there 
is  appreciable  loss  in  the  small  tunnel  dimensions  allowed  in 
HFSS.  As  we  can  see  from  Fig.  9(b),  there  is  strong  agreement 
in  the  axial  field  intensity  along  the  main  and  branch  tunnel 
axes  between  the  ADI-PE  and  HFSS.  The  figure  also  shows 
there  is  about  a  10  dB  drop  when  going  from  the  main  to  the 
branch  tunnel  (at  the  point  marked  C  in  Fig.  9(b)).  Although 
the  ADI-PE  is  used  to  simulate  a  tunnel  with  a  relatively 
small  electrical  cross-section  (2.7  Ax  3  A),  it  is  capable  of 
handling  larger  tunnels  at  higher  frequencies  without  running 
into  memory  problems  on  an  average  (2  GB  RAM)  PC  [1]. 

IV.  Rough-Walled  Tunnels 

A.  Rough-Walled  Tunnel  Model 

So  far  we  considered  only  tunnels  with  smooth  walls,  but  in 
this  section  we  investigate  the  effects  of  surface  roughness.  Sur¬ 
face  roughness  is  the  local  variation  of  the  tunnel  wall  relative  to 
a  mean  surface  level  [4],  [5].  In  this  study  we  consider  random 
surface  deviations  in  an  otherwise  smooth  wall.  For  the  purpose 
of  numerical  computations  we  assume  the  random  deviations  to 
be  Gauassian  distributed.  A  Gaussian  distribution  of  the  surface 
level  can  be  characterized  by  a  root-mean- square  height  devi¬ 
ation  ah  and  correlation  length,  Z[4],  [5].  Smooth  tunnels  have 
a  typical  RMS  height  deviation  of  0.01  m  and  rough  surfaces, 
such  as  those  found  in  coal  mine  tunnels,  have  a  RMS  height 
deviation  of  0.1  m  [4].  The  excess  loss  of  the  Ex  field  due  to 
roughness  in  a  rectangular  tunnel  is  given  by  (18) [4] 


Aough  =  4.3437T2(j|A2  ^2  +  z  (dB)  (18) 

where  di  and  d 2  are  the  tunnel  dimensions  in  the  x  and  y  axes, 
respectively,  and  2  is  the  range.  The  excess  loss  is  derived  by 
treating  the  rough  surface  as  a  random  process  and  from  ray 
tracing  techniques,  as  outlined  in  [10]  and  [4].  Equation  (18)  is 


shown  in  [4]  to  agree  with  experimental  data  taken  for  coal  mine 
tunnels  for  frequencies  ranging  from  200  -  1000  MHz. 

In  our  model,  we  treat  small  scale  roughness  by  replacing 
the  rough  surface  with  a  flat  impedance  surface  that  produces 
an  equivalent  specular  reflection  coefficient.  The  equivalent 
impedance  for  horizontal  and  vertical  polarization  is  shown  in 
(19)  and  (20),  [5] 


7  V  _ 
L eq  ~ 


zs  -  j(ko(?h)22j^i, 

Zs  +  (k0<Jh)2  sin  41, 

Zs  +  ( k0ah)2T () 

Z.-j(ko<rh)*>G0Z, 

Zs  +  (k0crh)2  sin3  1J', 


kQl  <  1 

kQl  >  1,^  >  yjp 

k°l  >  1,^  <  -^= 

(19) 

kj  <  1 
kj  > 

(20) 


where  is  the  surface  impedance  of  the  smooth  wall,  T(-)  is 
the  Gamma  function,  and  is  the  grazing  angle.  When  taking 
into  account  the  effects  of  roughness,  we  substitute  the  surface 
impedance  in  (3)  with  the  equivalent  impedance  of  either  (19) 
or  (20). 

Due  to  the  PE  angle  limitation  of  ±15°,  the  maximum  slope 
angle  0  of  the  rough  surface  and  the  grazing  angle  must  satisfy 
the  following  relationship 


20  +  ^  <  15°  (21) 

where  is  defined  as  shown  in  Fig.  10.  As  we  can  see  from 
Fig.  10,  the  angle  of  the  specular  reflection  of  the  incident  ray, 
denoted  by  £,  depends  on  the  height  deviation  of  the  roughness. 
The  roughness  angle  is  related  to  the  RMS  height  and  correla¬ 
tion  length  by 


tan  0  = 


0h 

l  * 


(22) 
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TABLE  III 

Analytical  and  Numerical  MAFs  for  the  Rectangular  Tunnel 


Rectangular  Tunnel  (dB/km) 

4.26  mx  2.10  m 

Analytical 

ADI-PE 

w/o  Roughness 

31.0 

29.8 

Roughness 

38.0 

36.1 

Excess  Loss 

7.0 

6.3 

TABLE  IV 

Analytical  and  Numerical  MAFs  for  the  Circular  Tunnel 


Circular  Tunnel  (dB/km) 
Radius  =  2.0  m 

Analytical 

ADI-PE 

w/o  Roughness 

1 1 .0 

10.1 

Roughness 

27.0 

26.1 

Excess  Loss 

16.0 

16.0 

B.  Comparison  of  Numerical  and  Analytical  Solutions 

We  consider  a  rectangular  4.26  m  x  2.10  m  tunnel  and  a  cir¬ 
cular  tunnel  with  radius  of  2  m.  The  fundamental  EHf  1  mode  is 
used  as  the  initial  field  of  the  rectangular  tunnel  and  the  funda¬ 
mental  TEoi  mode  generated  by  a  loop  ring  excitation  is  used  as 
the  initial  field  for  the  circular  tunnel.  Both  tunnels  operate  at  a 
frequency  of  1  GHz  and  the  dielectric  constant  and  conductivity 
of  the  tunnel  walls  are  taken  to  be  12  and  0.02  S/m,  respectively. 

Tables  III  and  IV  summarize  the  mode  attenuation  factors 
(MAFs),  or  the  loss  in  dB/km,  of  the  smooth  and  rough  tun¬ 
nels  with  rectangular  and  circular  cross-sections,  respectively. 
In  real  tunnels,  as  in  our  simulations,  the  lowest  order  mode  will 
determine  the  MAF  over  a  long  distance.  We  used  (18)  as  our 
analytical  loss  factor  for  both  the  rectangular  and  circular  tunnel. 
As  we  can  see  from  (18),  the  loss  due  to  roughness  is  a  function 
of  wavelength.  To  notice  an  appreciable  loss  at  1  GHz,  we  as¬ 
sume  the  walls  are  as  rough  as  cave  walls.  Therefore,  the  RMS 
height  is  0.1  m  (0.33  A)  and  correlation  length  is  2.5  m  (8.33  A) 
for  both  tunnels. 

The  grazing  angle  of  the  fundamental  mode  of  the  rectangular 
and  circular  tunnels  is  computed  using  the  analytical  expres¬ 
sions  for  the  propagation  constant,  (3  (where  e~^z  represents 
propagation  in  the  positive  z  direction),  outlined  in  [4]  and  [3]. 
The  grazing  angle  is  obtained  from  plane  wave  theory  by  [12] 


cos  A/  = 


Re(/3) 


(23) 


and  is  found  to  be  4.56°  and  5.25°  for  the  rectangular  and  cir¬ 
cular  tunnel,  respectively.  The  roughness  angle  is  found  to  be 
2.29°  from  (22),  and  (21)  is  satisfied.  The  ADI  is  simulated 
using  discretizations  of  Ax  =  0.284  A,  Ay  =  0.14  A  and 


Fig.  1 1 .  The  field  intensity  for  the  straight  tunnel  without  roughness  (solid),  for 
the  straight-curved-straight  case  (dashed),  and  for  the  straight-curved-straight 
case  with  roughness  (dot-line). 


A z  =  4  A  for  the  rectangular  tunnel  and  Ax  =  Ay  =  0.44  A 
and  A z  =  1.67  A  for  the  circular  tunnel.  As  we  can  see  from 
Table  III,  the  excess  loss  due  to  roughness  for  the  rectangular 
tunnel  is  about  7  dB  when  using  either  (18)  and  6.3  dB  when 
using  ADI-PE.  Simulations  using  a  unit  strength  Gaussian  initial 
field  (where  multiple  modes  are  allowed)  with  ax  =  ay  =  2.5  A 
and  ax  —  Gy  —  1  A  show  much  less  than  1  %  difference  in  MAF. 
Table  IV  shows  the  same  close  agreement  in  numerical  and  the¬ 
oretical  excess  loss  for  the  tunnel  with  circular  cross-section. 
In  this  case,  the  excess  loss  due  to  roughness  is  about  16  dB. 
The  accuracy  of  the  results  suggests  that  the  equivalent  surface 
impedance,  along  with  ADI-PE,  is  an  adequate  model  for  de¬ 
termining  loss  due  to  surface  roughness.  For  completeness,  we 
look  at  the  case  of  the  curved  tunnel  with  roughness.  Fig.  1 1 
shows  the  Ex  field  intensity  for  the  rectangular  tunnel  for  the 
straight  smooth  case  (solid),  the  curved  smooth  case  (dashed) 
and  the  curved  case  with  roughness  (dot-line)  with  a  radius 
of  curvature  of  1 800  m.  The  curved  section  is  in  between  two 
straight  sections  marked  off  by  two  vertical  lines.  The  MAF  of 
the  curved  section  with  roughness  is  about  37.7  dB/km.  Tunnels 
with  this  type  of  configuration  represents  a  topic  of  possible  fu¬ 
ture  work. 

C.  Comparison  of  the  ADI-PE  to  Experimental  Data 

In  this  section,  we  examine  a  real  underground  tunnel  with 
rough  walls  and  a  curving  axis.  The  tunnel  is  approximately 
3  m  x  3  m  and  the  walls  are  characterized  by  er  =  5.9,  a  = 
0.21  S/m  and  with  side  wall  roughness  of  ah  =  9.54  cm 
[11].  The  transmitter  is  placed  at  the  entrance  1.71  m  above  the 
ground  and  the  receivers  are  2.73  m  above  the  ground  and  op¬ 
erate  at  a  frequency  of  2.4  GHz.  The  data  was  recorded  at  every 
meter  along  the  length  of  the  tunnel.  The  curvature  of  the  tunnel 
axis  is  shown  in  the  boxed  region  in  Fig.  12.  The  experimental 
data,  shown  as  the  solid  line,  was  provided  by  the  research  team 
of  Moutairou  et  al.  [11]. 

The  numerical  simulations  were  done  by  computing  the  radii 
of  curvature  as  a  function  of  the  tunnel  axis  and  applying  rough¬ 
ness  for  the  side  walls.  A  unit  strength  Gaussian  wither^  —  ay  — 
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Fig.  12.  The  experimental  axial  field  intensity  of  the  rough-curved  tunnel 
(solid),  the  ADI-PE  simulations  (dashed).  (Boxed)  The  axial  geometry  of  the 
rough-curved  tunnel. 


2.5  A  is  used  as  the  initial  field.  The  correlation  length  was  se¬ 
lected  to  be  0.5  m  in  order  to  avoid  a  steep  slope  in  the  wall 
irregularities.  The  discretizations  are  taken  to  be  Ax  =  0.6  A. 
Ay  =  0.69  A.  The  equivalent  impedance  represents  averaged 
values  within  one  correlation  length,  so  As  is  chosen  to  be  2 
A  (0.25  m)  (not  much  smaller  than  the  correlation  length).  The 
ADI-PE  results,  shown  as  the  dashed  line,  are  shown  in  Fig.  12; 
In  the  simulations,  a  long  straight  section  (85  m)  precedes  the 
tunnel  in  order  to  minimize  the  presence  of  higher  order  modes 
in  the  ADI-PE  model.  A  sliding  average  with  a  window  of  1 
meter  was  used  to  match  the  resolution  of  the  experimental  data. 
As  in  previous  cases,  the  overall  trend  is  captured  after  the  near 
zone  region. 

V.  Conclusion 

The  ADI-PE  method  is  used  to  investigate  radiowave  prop¬ 
agation  in  three  specialized  tunnel  environments;  curving  tun¬ 
nels.  branch  tunnels  and  rough-walled  tunnels.  We  compared 
ADI-PE  numerical  examples  to  analytical  and  experimental  data 
or  to  numerical  data  provided  by  HFSS.  For  the  curved  tunnel, 
there  is  good  agreement  between  the  ADI-PE  and  known  analyt¬ 
ical  solutions  for  a  wide  range  of  curv  ature  radii.  For  branch  tun¬ 
nels,  even  at  the  PE  limit,  there  is  good  agreement  between  the 
ADI-PE  and  with  commerical  simulation  codes  such  as  HFSS. 
Also,  the  excess  loss  created  by  rough  walls  is  accurately  mod¬ 
eled  using  equivalent  impedances.  The  ADI-PE  method  com¬ 
pares  well  with  known  theoretical  loss  factors  for  tunnels  with 
either  rectangular  or  circular  cross-sections. 

For  each  case,  the  ADI-PE  method  accurately  models  propa¬ 
gation  in  the  tunnel  environment,  but  only  when  there  are  lim¬ 
itations  in  the  tunnel  geometry.  For  curved  tunnels,  the  curves 
must  be  smooth  and  propagation  step  sizes  of  2  -  3  A  is  ad¬ 
equate  for  accurate  results.  For  the  branch  tunnels,  the  branch 
angle  must  be  less  than  15°  to  ensure  accurate  results.  Within 
the  tunnel  junction,  where  the  staircase  approximation  is  used, 
propagation  step  sizes  of  less  than  1  A  is  recommended  for  ac¬ 
curate  results.  In  rough  tunnels,  the  maxirftum  RMS  height  de¬ 
viations  and  correlation  lengths  of  the  rough  walls  ire  limited. 


The  limitations  are  dependent  on  the  grazing  angle  as  well  as 
on  the  ratio  of  the  RMS  height  to  the  correlation  length.  In  this 
case,  simulation  data  suggests  propagation  step  sizes  of  about 
1 .5  -  4  A.  As  shown  in  previous  work,  the  discretizations  along 
tlie  x  and  y  axes  must  be  less  than  1 .9  A  to  satisfy  Nyquisl’s  the¬ 
orem  [1 1.  However,  for  the  cases  considered  here,  good  results 
were  obtained  for  discretizations  less  than  1  A  for  the  curved  and 
rough-walled  tunnels  and  less  than  0.5  A  for  the  branch  tunnel. 
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